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ABSTRACT 



An algorithm Is developed for estimating the state of a linear dy- 
namic system excited by a random sequence. The input data are noisy 
observations which are nonlinear functions of the state. The estimates 
are best in the sense of least squared residuals. A significant problem 
in radar tracking is investigated and the effectiveness of the algorithm 
verified. 
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1 . 



Introduction. 



The problem to be considered is that of state estimation where the 
observations are nonlinear functions of the state of the system corrupted 
by additive white noies. The equations of motion of the state are lin- 

V 

ear in the state and the excitation which includes random components. A 
significant constraint on the solution is that the conceptual solution 
to this problem must lead to an explicit procedure which can be realized 
on a digital computer and the scheme must produce estimates as the obser- 
vations are received. 

The theory for the case where system dynamics and observation func- 
tions are linear is very highly developed. However, when nonlinearities 
are introduced there are virtually no completely satisfactory solutions. 

Two methods for handling the nonlinear problem have been introduced. 
The first entails a linearization about a nominal trajectory in state 
space. Its success depends upon the accuracy of the nominal trajectory. 
This technique has little hope of success in a situation where there is 
almost no prior information about the trajectory. The target acquisition 
problem is an example of such a situation. 

The second method is more of theoretical interest than as a candi- 
date for a computational procedure. In this development the viewpoint is 
taken that the output of a state estimator should be the conditional pro- 
bability density function, conditioned upon all past data. Computational 
difficulties arise from an effort to compute a complete function over the 
entire state space as compared to a more conventional estimator which 
selects a single point in the state space as the most likely state. 

This study was particularly motivated by the difficulties encountered 
in filtering radar returns from an airborne target. The target dynamics 
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are presumably descrlbable by a linear dynamic system where the elements 
of the state vector are the position and velocity of the target in car- 
tesian coordinates measured from the radar. On the other hand, data avail- 
able to a filter from the radar will usually be in spherical coordinates. 
Thus there exists a known, nonlinear relationship or transformation be- 
tween the state of the system and observations. The results of a series 
of experiments using a filter based upon a simple linearization procedure 
are reported in Demetry and Hudson [4]. The operation of the filter was 
unsatisfactory under realistic conditions of initial uncertainty about 
target position, i.e., where to evaluate the partials involved in the 
linearization procedure. 

The present work fills a gap in the field of nonlinear estimation for 
problem in which there is little prior Information and a computationally 
feasible estimation procedure is required. 

The problem is discussed and precisely formulated in Section 2. The 
original filtering problem is replaced by an associated minimization prob- 
lem. The work of Kalman is reviewed since it is known that the Kalman 
filter is the solution to the associated minimization problem when the 
observations are linear functions of the states. A special single-stage 
form of the problem is considered. It is shown that the nonlinear obser- 
vation problem can be approximated by a linear problem whose solution is 
known from Kalman's work. The resulting solution is then used to generate 
a new approximation to the original problem. Each iteration of the above 
process reduces the function to be minimized so long as the gradient is 
non-zero. It is shown that the whole class of problems originally con- 
sidered can be cast in the special single-stage form. 

The problem may be said to be solved at this point; however, for real- 
time calculations there must be procedures for controlling the number 
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of iterations. Iteration control procedures are discussed along with 
heuristic means for choosing the control parameters. Even with the num- 
ber of iterations kept to a finite numberi the computing requirements 
would increase indefinitely since each new set of data generates a new 

3 '• 

minimization problem over a larger number of variables. The control of 
this effect is discussed by introducing the concept of noise generated 
by the nonlinearities. The overall algorithm is discussed with respect 
to implementation on a digital computer. 

The method was used on a realistic radar tracking problem. The 
models for the target and the radar are discussed. The results of the 
study indicate that the method produces reasonable estimates of the states 
of the system. 
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2 . 



Detailed Statement of the Problem. 



The specification of the problem may be viewed as consisting of four 
parts. Three of these are different types of information about the system 
under observation and the last is an implicit statement of how the avail* 
able information should be combined to form an estimate of the state of 
the system. 

The first type of information about the system is expressed by form- 
ing a dynamic model of the system. In this way the relation between 
states at different times is made explicit. It is only through the dy- 
namic relation of the states that observations taken at diverse times 
have any relation to one another. In general usage the term filter is 
associated with a sequence of observations which are related to one 
another (correlated). The states and the dynamic model embody this cor- 
relation. 

The second type of Information is the relationship of the observa- 
tion at a given time to the state at the same time. 

The third type of Information is the a priori knowledge about the 
state of the system. 

Finally the best estimate of the state is defined. In general, none 
of the information about the state is definitive; it is all subject to 
some uncertainty and except for this uncertainty, it would be contra- 
dictory. The best estimate is defined to effect a certain compromise 
among all of the available information. 

Since the resulting definition is implicit, computational difficul- 
ties occur. The resolution of these difficulties results in an explicit 
procedure for producing an estimate which is almost best within a reason- 
able computing time. Such an explicit procedure, or algorithm, will be 
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called a filter. 



The Dynamic Model 

The time evolution of the states of the system which is being observed 
are assumed to be adequately described by a difference equation which is 
linear in both the state and the excitation. The excitation is assumed to 
be a white random time sequence which has known first and second moments 
and is independent of the state. If the mean of the random excitation is 
not zero, then the random signal can be decomposed into a deterministic 
component and a zero*mean random component. The development will assume 
that there is no deterministic component (or non-zero mean). A short com- 
ment will be made at the appropriate point outlining the required changes 
for the case of deterministic inputs. These notions are concisely stated 
in (1) through (4). 



- 4 X(k) + r w(k) 


(1) 


Eru)(k)] ^ 0 


(2) 


Erui(k) . (J IjJ 


(3) 


E[uj(k) x(j)^] = 0 for k^ 


(4) 



where x is the state vector of n components, 

r is a known n x r input distribution matrix, 
uj is a random excitation of r components, 
is a known r x n state transition matrix, 

Ef ] is the expectation operation, and 
T denotes the transpose. 

Equation (1) expresses the linearity of the state dynamics while (2), (3), 
and (4) express the qualities of zero mean, whiteness, and independence 
respectively about the excitation. The assumption that the covariance 
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of Ui)(k) Is the identity matrix involves no loss of generality as long as 
r is appropriately chosen. Consider the two random excitations I’d) and 
r'ti)' where 

E[uj] » 0 
Eru)'] = 0 
Efuj U)*^] = I 
E[u)' = Q 

By comparing first and second moments, the two random excitations are 
X T 

equivalent if F F = F'OF'^* Q is a covariance matrix and thus is sym- 
metric and positive semi-definite. This implies that a decomposition can 

T 

be found such that B B = Q, and F •• F'B. 

The Observations 

The data available to the filter are nonlinear functions of the 
state corrupted by additive white noise. The functional relation is 
assumed to be twice differentiable in the state. The corrupting noise 
is assumed to have zero mean and known variance. The noise is assumed 
to be independent of the states and the excitation. These notions are 
concisely expressed as 



ir(k) = h(x(k)) + U(k) 


(5) 


Eru(k)] m 0 


(6) 


EfPOc) III IJJ 


(7) 


E[u(k) = 0 


(8) 


Eru(k) = 0 


(9) 



where ^r(k) is an m vector of observations at time k, 
y(k) is an m vector of noise at time k, 
h( ) is an m vector of nonlinear functions of the states 
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M Is the m X m covariance matrix of the measurement noise* 



A Priori Information 

In view of (1) any information relative to Je(k) must also be consid- 
ered when estimating JC(k4“l) . In conventional linear sequential stage-by- 
stage estimation one can consider two distinct phases. The first is to 
bring forward all information from the past observations in the form of 
an a priori estimate using (!)• The second phase is then to adjust this 
a priori estimate in view of the actual observations of the state. This 
is repeated from stage to stage. Clearly this process must start with a 
given a priori estimate for the first state. Sometimes this a priori 
estimate serves merely as a mathematical convenience for starting the fil- 
ter [ 9 J. In other cases an appropriate a priori estimate is really avail- 
able. In any case an a priori estimate takes the form of an estimate of 
the initial state coupled with a measure of the accuracy of this estimate, 
the covariance of the error, defined as 

E[Cr(l)-a;(l/0)) (at(1/0))'^] s P(1/0) (10) 

E[jr(l)-Jf(l/0)] » 0 (11) 

where x(l/0) is the a priori estimate of JC(1), and 

P(l/0) is the covariance of the a priori estimate. 

The double index argument will be used throughout to indicate an estimate 
of the state associated with the first index based on observations up to 
and including those associated with the second index. 

Definition of Best Estimate 

The best estimate of sequence of states x(l) through :r(k) will be 
called jc(l/k) through x(k/k) and is defined as that sequence which mini- 
mizes the scalar quantity 
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CCx(l/k), jc(2/k) x(k/k)] s || x(l/k) - jf(l/0) 

Jc 

+ ^ II ;e(i/k) - ♦ard-l/k) ||J 

i*=2 ^ 

_k 

+ ^ II <?(i) - hC3ff(i/k))||^ (12) 

1=1 ^ 

where Che norm notation is introduced for compactness. By definition 
0 x 

||b|K is equivalent to b Ab and the result is a scalar which is a quad* 

A 

racic function of Che elements of b. 

A best estimate defined in this manner might be called a weighted- 
least-squares estimate generalized to include the case where the quantity 
Co be estimated is changing somewhat randomly in time. 

The three different types of terms correspond to (10), (1), and (3) 
respectively. The terms inside Che vertical bars are called residuals. 
The residuals are the difference between the expected value of a function 
of Che true states and chat same function evaluated at Che estimate of 
Che state. 

This definition of Che best estimate can be interpreted in several 
ways which will be developed below. These interpretations cannot in any 
sense prove that estimates defined in this way are best. The most that 
can be hoped for is that the interpretations offered will enhance che- 
reasonableness of the resulting estimate. It must be realized that Che 
best estimate is only what it is defined to be, which, in Che final anal- 
ysis is certainly somewhat arbitrary. 

First, if all of the random sequences are assumed to be sequences of 
Gaussian (or normal) random variables, then it is possible Co compute the 
probability of the observed data for a given sequence of states. This 
probability is viewed as a function of Che true states. That sequence of 
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the states which maximizes the probability is taken as an estimate of 
the true states. The probability is commonly called the likelihood and 
can be expressed as 

L(Jf(l), 0f(2), Jc(k)) ■ K exp [- C(jf(l), J:(2) , Jf(k))] 

where K is a constant which is independent of the states if the weighing 
matrices are chosen as 

W^-Pd/O)'^ (13) 

W^-Cri^)'^ (14) 

(15) 

The covariances are assumed to be nonsingular for simplicity. For 
a detailed discussion of the case of singluar covariance matrix see 
Appendix I. It is clear that to maximize L one must minimize C. Thus 
for the case of Gaussian random variables the best estimate will be the 
so called maximum likelihood estimate. 

A second point of view is that it would be desirable to find an 
estimate which resulted in zero residuals; requiring zero residuals » how- 
ever, would imply a larger number of constraints than there are adjustable 
parameters (estimates). This suggests that the best estimate would be 
some sort of compromise where all of the residuals are small. Such a com- 
promise is effected by setting up a weighted sum of squares of the re- 
siduals, C, as a function of the estimates, and selecting (or defining 
in this case) that set of estimates which minimize C as the best esti- 
mate. 

In general, the best estimate will depend upon the weighting chosen 
for each residual. In order to determine the appropriate weighting mat- 
rices it is helpful to consider under what circumstances equal weighting 
would be appropriate. A heuristically reasonable answer might be to 
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weight equally when each of the random components have the same variance. 
If each of the sets of equations are multiplied by appropriate matrices 
new random variables can be defined so that each has unit variance. The 
residuals from these new equations are computed and their squares all 
weighted equally. In this form the residuals have an intuitively rea- 
sonable weighting. But it can be shown that this is equivalent to 
weighting the original residuals with the Inverse of the corresponding 
covariance. 

A very simple example should clarify the argument. Consider the 
problem of estimating x from two observations 

Z ” hj^Of) 

where 

Ery2] = o» 

EfWi] = 

EfU^] = 

Dividing the first equation by the second by a£ yields 

‘S’l/t^l ■ h^(jf)/a^ + 

and 2 “ •'2 ^2 

where u 2 and U 2 are unit variance random variables. So the sum of squared 
residuals is formed 

c(x) * + (*2/(^2 “ 

but this is equivalent to 

COf) - (* (»2 - h20c))^/a2 . 
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Therefore, if it is reasonable to weight equally residuals associated 
with random variables of equal variance then it follows for unequal vari- 
ances the weighting should be in inverse proportion to the variance. 

Another point of view is that the weighting matrices Wj^, W 2 and 
are chosen directly without recourse to any assumed random variables. 

There is, of course, an equivalent problem cast in terms of random vari- 
ables. It is the author's view that it is easier to assess the magnitude 
of the variances of the random variables than to choose an appropriate 
set of weighting matrices directly. 

A general model of the underlying physical process which generates 
the data has been presented. For any real physical situation the para- 
meters r, x(l/0) and P(l/0) will be numerical quantites and the non- 
linear functions h(x) will be a vector of explicit functions. This is the 
type of information which a filter designer must have before the filter 
can be constructed. For a given model there are many possible filters 
which might be considered; in each case presumably the output of the fil- 
ter would be best in some sense, often unspecified. The filter under 
consideration in this work is defined by the sense in which the estimates 
are best, i.e., the filter is a least-squared-residuals filter. Thus, it 
should be noted that the filtering problem has been transformed into a 
sequence of minimization problems. It will turn out that the solution to 
each minimization problem is itself a sequence of solutions to a much 
simpler minimization problem, namely the problem with linear observation 
functions. The solution to this simpler problem is well known and is 
discussed in the next section along with other filter techniques where 
the model is similar to the one described in this section. 
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3. 



Prior Work. 



The field of state estimation has received a great deal of interest 
recently, especially since the work of R. E. Kalman [7] and R. E. Kalman 
and R. Bucy C8] in linear estimation theory. For a comprehensive survey 
of the general field of estimation Deutsch f6] is suggested. Lee [9] 
presents a fine treatment of the theory, particularly with respect to the 
relationship between control and estimation theory. Cox f2, 3] is a 
specialized review of the efforts in the area of nonlinear estimation of 
which this work is a special case. 

The structure of the problem and the results obtained by Kalman [7] 
along with several extensions, modifications and alternate interpretations 
are discussed in some detail. This discussion provides a convenient refer- 
ence for comparison of the results of this thesis as well as an opportunity 
to establish certain known results which will be needed in the development 
of the method that follows. The method of trajectory linearization is dis- 
cussed and it is shown how a natural extension leads to the method used 
here. 



Kalman Filter 

Kalman [7] has solved a special case of the problem under consider- 
ation where the measurements are linear functions of the state. 

Symbolically 

2(k) - F JC(k) + U(k) (16) 

replaces (5). 

Two cases are considered. In the first, all of the random sequences 
are assumed to be sequences of Gaussian random variables. With this as- 
sumption, the estimate is shown to be optimum in the sense that any linear 
function of the estimate is the minimum variance estimate of the same 
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linear function of the true state. The estimate turns out to be a linear 



function of the observations. 

In the second case, there are no assumptions about the form of the 
probability density functions of the random variables, but the estimate 
is assumed to be a linear function of the observations. The optimum 
estimate is defined in the same way. 

In either case the method of computing the optimum estimate (the 
filter) is the same. The sequence of operations can be envisioned as 
consisting of two steps. The first will be called the prediction equation* 

jc(k^-l/k) » $jc(k/k) (17) 

The double argument notation will always indicate an estimate. The left 
side should be read; the estimate of the state at time given data up 
to time k. The second step will be called the adjustment for newly re- 
ceived data. 

x(k+l/kfl) = xOcfl/k) + G(k) [^(kfl) - H A:(kfl/k)] (18) 

where [^(kfl) - ^ j;(kH/k)] is the error in the predicted observations, 
and(7(k) is a matrix of adjustment coefficients. The matrix C?(k) reflects 
the relative confidence one should have in the observed data as compared 
to the predicted estimate. This is discussed in 



G{k) = P(kfl/k) r^P(kfl/k) (19) 

Where P(k+l/k) is the covariance of the estimates defined as follows: 

P(k+l/k) = E[Cc(k+l/k) - .r(kfl))Cx(k+l/k) - x(k^-l)^] . (20) 

This formulation then requires that one must keep track of the co- 
variance of the estimates. This is also done in two steps. 



p(id-i/k) = ♦ p(k/k) 4'"’^ + r 

p(k+i/kfi) = 

P(kf l/k) - p(kf l/k) IH P(kf l/k) + p]"^ H P(kf 1/k) 



( 21 ) 

( 22 ) 
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Rauch [13] has derived these same equations based on the Gaussian 
assumption and shown the resulting estimate is the conditional mean and 
the maximum likelihood estimate* 

Lee [9] has shown that the same equations yield a weighted least 
squared residual estimate where the weightings are the inverses of the 
covariance matrices. This also follows from the derivation of the maxi- 
mum likelihood estimate based upon the assumption that the random se- 
quences are Gaussian. 



Trajectory Linearization 

The trajectory linearization method has been used to solve nonlinear 
filtering problems such as orbit determination for artificial satellites 
[10], [11]. It is assumed from physical considerations that the evolution 
of the system state satisfies a general difference equation of the form 

= fCc(k), Uj(k), k) (23) 

and that observations are available in the form 

^(k) = hOc(k), U(k), k) (24) 

where u}(k) and l^(k) are random vectors. 

It is further assumed that there is a known nominal trajectory which 
is a sequence of states jc°(k) such that 

jc°0cfl) = fO°(k), 0, k) . ^'25) 

It is desired to find the best estimate of the deviation of the true 
state from the nominal state. 

Defining 

^(kfl) a xik+1) - jc°(k+l) (26) 

and substituting (23) and (25) into (26) results in 

y(k+l) = f(x(k), u,(k), k) - fCc°(k), 0. k) . (27) 

This relation is now approximated by a first order Taylor series about 
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the point x(k) = ;c°(k) and (u(k) ■ 0. The two partial derivatives are 
given appropriate symbols. 

4 (k) . f 0 c°(k), 0 , k) (28) 

JC 

T (k) s 0 , k) (29) 

The first order Taylor series expansion results in 

^(kfl) ^ fCx°(k), 0 , k) + 4(k) U(k) - x°(k)] 

+ r(k) uj(k) - fCc°(k), 0, k) . (30) 

Substitution of (26) in (30) results in 

« ^(k) j/(k) + r(k) u)(k) . (31) 

From the nominal trajectory it is possible to construct a nominal 
set of observations 

2 °(k) = hCc°(k). 0, k) . (32) 

Consider the deviations of the observations about the nominal obser- 
vations . 

U(k) s ^(k) - ^°(k) (33) 

Substituting (24) and (32) in (33) yields 

u(k) = hCr(k), y(k), k) - hCr°(k), 0, k) (34) 

The relation is approximated by a first order Taylor series about 
the point ;c(k) = o:°(k) and V(k) = 0. The two partial derivatives are 



given appropriate symbols 

^(k) 5 h^Cr°(k), 9 , k) (35) 

S(k) B hj^(x°(k), 0, k) (36) 

u(k)~ h(je°(k), 0 , k) + ¥(k) roc(k) - :t°(k)] 

+ S(k) U(k) - hCr°(k), 0, k) (37) 

Substitution of (26) in (37) yields 

u(k) = ^(k) y(k) + S(k) y(k) (38) 
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Although it was not noted in the description of the Kalman filter, 
it is true that the equations remain valid if any or all of the matrices 
H, T, B, are known functions of time. These filter equations are thus 
directly applicable to (31) and (38). 

The purpose of the process of trajectory linearization is to generate 
(31) and (38). It is then noted that with respect to the states j/(k) and 
the observations u(k) the model is in the form of a linear dynamic system 
and linear observations. The Kalman filter is then applied directly as 
though (31) and (38) were equalities. 

Nonlinear Noise 

A question naturally arises concerning the adequacy of the first ord- 
er approximation in developing (31) and (38). The heart of this problem 
is investigated by Denham and Pines f 5] through the use of a very simpli- 
fied model and a number of Monte Carlo studies. They reach the conclusion 
that the difficulties are of an indirect nature. The first estimates are 
about as good as might be expected. In processing subsequent data, how- 
ever, trouble develops because the assumed quality of the first estimate 
is too great, which means that the next data get weighted too lightly. 

This effect can be best seen by reconsidering (31). In order to make 
this expression into an equality, all of the higher order terms must he 
added to the right side of the expression. These additional terras should 
be considered as part of the observation noise. It is the failure to 
account for this nonlinear noise in the Kalman filter that causes the dis- 
crepancy between the covariance of the estimates as computed in the filter 
and the true average squared estimation errors. When the calculated co- 
variance overstates the quality of a given estimate, a subsequent obser- 
vation will certainly be combined with the current estimate in a non- 
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optimum manner. 

Denham and Pines point out that when the order of magnitude of the 
expected value of the neglected terms is of the order of the natural 
measurement noise one cannot expect the linearized filter to work proper- 
ly. The expression for the **nonlinear noise" involves the difference 
between the true state and the point about which the linearization takes 

place. If there were some way to reduce this difference then the nonlinear 

noise would be reduced correspondingly. These authors attribute to John 
Breakwell an iterative procedure for accomplishing this. The procedure 
is to linearize and filter, then relinearize at the new estimate and fil- 
ter again. This cycle is repeated until the output of the filter is the 
same as the point at which the linearization takes place. A Monte Carlo 

study using this iterative technique revealed that the computed covariances 

quite accurately reflected the quality of the estimates. 

It will be noted that the method used to solve the least-squared -re- 
sidual problem as developed in the next section is exactly the iterative 
method suggested by Breakwell. 
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4. Development of the Solution Algorithm. 

The development will proceed in two phases, the first being a con- 
ceptual means of finding the absolute best estimate, the second being the 
development of a series of compromises required for computational feasi- 
bility . 

The minimization of (12) will have to be carried out in an iterative 
fashion since the simple process of differentiating and setting to zero 
does not lead to an explicit formula for the state estimates as it would 
if the observation functions were linear. The iterative procedure is 
based upon a linearization of the observation functions. The linearized 
observations are then in the form (16) and minimization is carried out 
using the Kalman filter equations. This produces a set of state estimates 
about which the nonlinear functions can be relinearized. This process is 
repeated until there is no further change in the state estimates. 

The Kalman filter in its normal form is not completely adequate since 
its output is the sequence of estimates jc(l/l) through jc(k/k) while the 
point about which it is desired to linearize is ,y(l/k) through x(k/k) . 
These latter estimates are called the smoothed estimates. There are for- 
mulas available, due to Rauch [13], for converting the output of the Kalman 
filter into smoothed estimates. There is, however, a more convenient way 
to get the same results in this case where all of the smoothed estimates 
are required. This involves converting from a multistage problem to a 
single-stage problem with a proportionately enlarged state space. The de- 
tails of this conversion will be discussed following the discussion of the 
iteration for the single-stage process. 
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5. The Single Stage Minimization Procedure. 

The equations related to the single-stage process are rewritten here 
using a simplified notation and are numbered using the same numbers as in 
their first appearance with an * added. 



2 ~ h(x) + V 


(5') 


E[U] = 0 


(6') 


E[uu’’^] = P 


(7') 


EfUjc^] > 0 


(S') 


Ef Oc-jCq) = P^ 


(10') 


= 0 


(11') 


1 - *JI p * 11^ - 

o 


(12') 



where 

Z is a vector of observations , 

U is the noise in these observations, 

R is covariance of the noise, 

is the a priori estimate of the true state x* 

is the new estimate of Xt and 

is the covariance of the a priori estimate. 

The pertinent equations from the solution to the linear problem are 
also rewritten here , 

z = H X ^ V (16') 

jCj = + G(z-Hx^) (18') 



G = 



( 19 *) 



The Don>slngularlty of P^ assumed in (12*) is fully discussed in 
Appendix 1. It is sufficient here to say that if P^ is singular the 



27 



problem can be reduced to a smaller state space where the associated 
is non-singular. It will always be assumed that ^ is non-singular, i.e., 
the assumption is made that there are no observations of unlimited ac- 
curacy. 

The iterative minimization procedure involves a linear (first-order) 
approximation to (5') about a point which will always be the best 
available estimate of x* This linearized approximation is then manipulat- 
ed to form a synthetic observation which has the form (16*). This syn- 
thetic observation is used in (18*) and a new approximation to the best 
estimate is obtained. Using this estimate as the point of linearization 
of (5*) the process is repeated. These steps are expressed symbolically 
as follows. 

2 = h(x[) + + u (39) 

where jfJ is the ith approximation to the best estimate x. 

2 - hOfJ) + jeJ (40) 

+ U (41) 

where is the so-called synthetic observation and 

= h^Ofj) (42) 

+ X^) (43) 

where 

(44) 

This completes the description of the pure minimization algorithm 
except for a specification of the initial point of linearization and the 
possibility of overshoot. Discussion on the initial point of linearization 
will follow after the discussion of conversion from multi-stage to single 
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stage. The possibility of an overshoot results from the fact that a 
simple first order approximation to the nonlinearity may not be accurate 
for the subsequent change in The overshoot protection scheme used 

is simply if 

a C(jrJ) 

then is replaced by 1/2 X^) * 

It will now be shown that each step in this process does result in 
a decrease in the cost, C. Since if the new cost is greater than the old 
cost the step size, > is reduced by a factor of 2 it only is 

necessary to show that for a small enough step size there will be a re* 
duction in C. The demonstration will be begun by showing that the change 
in estimate is related to the negative gradient of the cost evaluated at 
by a positive definite matrix. 

gCxJ) = (x[) (45) 

g(xj) = P'J' (46) 

The change in the iterate is 



. _ i+1 i 

d = 


(47) 


d = jf - 

O 1 o 


(48) 


d a - jcJ + (/*‘[<ar-hCx^) + 


(49) 


d = ri-(?V]0e^-;fJ) + (7^r^-h(xJ)] 


(50) 


show that 




d = P^P^I g(x[) 


(51) 



it must be shown that 
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[I - ffV] = CV^/" 



(52) 



and that 

+ P]"^ pVq] . (53) 

The first of these is obvious after substituting for , 

In the latter P is factored from the left side to obtain 
o 

+ P]'^ P’^. 

and then + P]”^ TpVqP^^ + P] is substituted for the I above. 

[pVqP^"^ + P]'^ \H^PjJ^'^ + P - pV/^^] P’^ 

After cancellation, the above is the expression for • 

Note the matrix + /?] ^ covariance of 

the updated estimate in the linear case so it will be given the special 
symbol 

P^ H ?o " P/^^ + P]’^ pVq (54) 



so that (51) can be written 



d = P^gOfj) (55) 

After applying the overshoot control the actual change of the iterate 
is in the direction of d but may have a smaller magnitude. Let D be the 
actual change so that 

D = qd (56) 

D= qP^gOj) (57) 

where q is some integral power of (1/2). 

Then to a first-order approximation the change in C, ^C, is given by 



AC» 2g^D 
AC= 2q g^'^d 



(58) 

(59) 
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(60) 



AC» 2q g^Pjg 



AC= 2q ||g|| 



2 



(61) 



Since all higher order terms in the expansion of AC involve higher powers 
of q it can be asserted that for some small enough q the first order term 
will dominate. Thus for some suitable q the change in C will be negative 
for any non-zero g if is a positive definite matrix. is known from 

the linear theory to be positive definite for any value of and any 
positive definite P. This also follows from the fact that 



a convenient matrix identity discussed in Tl], and the fact that the in- 
verse of a positive definite matrix is positive definite. 

It has been shown that the special single-stage problem can be solved 
by solving a sequence of simple problems which approximate more and more 
closely the real problem. Each iteration was shown to reduce the cost 
unless the gradient of the cost was zero. 

In the next section it will be shown that the general problem con- 
sidered in this method can be recast in the form of a single-stage prob- 
lem. 




(62) 
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6. Multi-stage Case Cast in Single stage Form. 

The process of converting a multi-stage estimation problem into a 
single stage problem is accomplished by defining the new state to be the 
juxtaposition of the states at the various stages. The process will be 
carried out in detail for a two- stage problem and then the process will 
be generalized by induction for a k-stage process. 

Let the new state be 



X = 



;e(i) 

mmmm 

X(2) 



(63) 



let the new estimate be 



X 



1 ■ 




9 



(64) 



let the a priori for the new state be 




and let the a priori covariance of the new state be 



P 



o 



>(1,1/0) I P(l,2/0) 
P^(l,2/0) 1 P(2,2/0) 



(65) 



( 66 ) 



Some of the elements in and p^ have not been previously defined. 

However, the notation used has already been defined, i.e., x(2/0) means 

the estimate of the second state given no data. The subraatrices in P ■ 

o 

are defined as follows: 

P(i,j/0) = E[Cr(i/0) - x(i»Ce(j/0) - .»(j))'‘’] (67) 



This is a natural extension of the double argument notation alreay defined 
which is necessary to accomodate consideration of the cross correlation 
between states at various times. Since these new elements are not given 
directly in the multi-stage model it will be necessary to fill in these 
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elements in order to define the single-stage model. 

First consider jc(2/0) . Based upon the requirement that (11) be sat- 
isfied it must be true that 

jc(2/0) == E[jr(2)] . (68) 

Substituting for jc(2) from (1) and taking advantage of (2), the fact that 
0 )( 1 ) has mean zero, yields 

x{2/0) = $E[o:(l)] . (69) 

Introducing (11) above yields 

x(2/0) = 4x(l/0) . (70) 

If there are deterministic inputs (or equivalently Er«»(k)] ^ 0) the 

appropriate modification is 

je(2/0) = ?';f(l/0) + r Eru)(l)] . ( 71 ) 

Now consider the various submatrices of p^, p(l,l/0) is already 

known in the multi-stage problem as P(l/0). p(l,2/0) is given by 

P(l,2/0) = EfOcd/O) - :e(l))0f(2/0) - jc(2))'^] . (72) 

Substituting (70) and (1) in the last part of (72) yields 

P(l,2/0) = Er0e(l/0) - je(l))(e;e (1/0) - - Toud))'^] . (73) 

T 

Taking advantage of the independence of (u(l) from (2) and factoring 
yields 

P(l,2/0) = E[0c(l/0) - jf(l))Ce(l/0) - :c(l))^]4>’^ (74) 

but this is just 

Pd, 2/0) = P(l/0)4>‘'^ . (75) 

By similar arguments it can be shown that 

P(2,l/0) = 4>P(l/0) (76) 

and 

p( 2 , 2 /o) = 4>Ki/o)<r’’^ + rr’’^ . (78) 

The remaining elements needed to complete the description of the single- 
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stage model have to do with the observation process 



Let the actual observations be 



■H- 



( 79 ) 



let the vector of observation functions be 

h(je(l))| 

- - - J ^ 

_hCr(2))| 



h(x) = 



(80) 



let the measurement noise be 



..M . 

L«(2)J 



(81) 



and let the measurement noise covariance be 






'r I 0 

l" 



(82) 



_0 I ^ 

where the off-diagonal elements of R are zero matrices in view of the 
whiteness of the observation noise. 

As the number of stages is increased, the dimension of the single- 
stage X Is also increased. The process of expanding the single-stage 
model of a k stage model to that of a (k + 1) -stage model is 
explained in detail only for Xt and The expansion process or aug- 

mentation for the remaining elements of the model is as simple as the aug- 
mentation of X will prove to be. The augmentation of x is shown explicitly 
as an example. 

In the augmentation process x goes from 






X = 



Jf(k) 



to 



'jc(l) 



JC(k) 

;r(k + 1) 
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The expansion of is only slightly more involved. For the (k 1)* 



stage case the a priori estimate is 

xH/0) 



*o = 



(83) 



jc(k/0) 

_x(k+ l/0)_ 

where + 1/0) = 4yc(k/0). 

The structure of the corresponding to the (k + l)-stage case is a 
(k + 1) by (k + 1) square matrix whose elements are submatrices. The 
upper left k by k part of this matrix is already known from the assoc- 
iated with the k-stage model. Thus to expand to the k + 1 case it is only 
necessary to fill in the lower border. The a priori covariance P^ has the 
form 



^o = 



P(l,l/0) 

• 


. . P(l,k/0) 

• 


1 

1 

1 


P(l,k+ 1/0) 
• 


• 

P(k*l/0) 


• 

• 

. . p(k,k/0) 


1 

1 

1 


• 

• 

p(k,k + 1/0) 


P(k + 1,1/0) . 


. . P(k + l,k/0) 


1 


P(k + 1, k + 1/0^ 



(84) 



The formulas below for the new border elements of P were derived in ex- 

o 

actly the same way the submatrices P(l,2/0), P(2,l/0) and P(2,2/0) were 
found in the case for k -h 1 =: 2. 



P(k + l,i/0) = <tP (k,i/0) for 1 ^ i ^ k (85) 

P(i,k + 1/0) s P*^(k + l,i/0) for 1 ^ i s k (86) 

P(k + l,k + 1/0) = (k,k/0) + r (87) 

Thus it is clear that the dynamic relations represented by (1) for the 
multi-stage problem are incorporated into the very special structure of 
the large covariance matrix P^ in the single-stage equivalent. 

After a given multi-stage problem has been cast in single-stage form. 
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the single-stage problem is solved using the methods previously described. 
The outcome of the single-stage solution is the best estimate, which 
must then be interpreted in terms of the original multi-stage problem. 

The best estimate, is the best estimate of all states given all data 
up to the present stage, i.e. 



Jf(l/k) 



jf (k/k) 



( 88 ) 



Finally, as each new single-stage solution process is started there 
must be an inital estimate of which is called jfJ . The first iterate 
for a (k -f l)-stage minimization will be based on the final estimate from 
the previous minimization over k stages. Explicitly, the first iterate is 
given by 

Jf(l/k) 



where 




jp(k/k) 

Jf(k + 1/kJ 
jf(k + 1/k) = <J>r(k/k) 



( 89 ) 



( 90 ) 



At this stage in the development of the filter, the problem has been 
solved but in a very impractical way. The filter has been broken into an 
unlimited sequence of minimization problems. Each minimization problem 
is solved through an as yet unlimited sequence of approximate solutions. 

In order to design a practical filter a realistic convergence test must 
be used to terminate the minimization process. Such a convergence test 
is discussed in the next section. Another difficulty arises in connection 
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with the sequence of minimization problems. As more and more data are 
considered, spanning a larger and larger collection of states, the size 
of the minimization problem increases. A practical means of limiting 
the size (dimensionality) of the minimization problem will be discussed 
in a subsequent section. 
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?• Criteria for Teminatlon of the Minimization Procedure. 

It seems to be a universal feature of any iterative numerical method 
that the termination decision is based upon **feel*^ or **rules of thumb**. 
Typically, a quantity is chosen as a measure of the convergence of the 
iterative process and this measure is compared against a standard or 
threshold. Fortunately, in this particular case it is possible to offer 
some insight into the choice of both the measure and the standard or 
threshold. This is true because the problem is basically a stochastic 
one* By analogy with certain special cases it is possible to approximate 
the probability density function of the measure of interest. 

The control law or algorithm for the control of the number of itera- 
tions is based upon the fact that the minimum cost, is Itself a 

random variable whose distribution function can be approximated by that 
of a chi square random variable. The minimum cost is exactly distributed 

as a chi square variable when the measurements are linearly related to 

the states and all of the random variables are Gaussian. The chi square 
distribution is characterized by a parameter called the number of degrees 
of freedom. For the least square problem the number of degrees of free- 
dom is the number of constraint equations less the number of parameters 
adjusted in the process of minimizing the sum of the squares. 

Using this assumed distribution function for the minimum value of C 

it is possible to evaluate numerically the probability of a minimim C being 

greater than some number C . Actually the question is reversed so that a 

C- is found such that the probability that a minimum C is greater than C. 
b b 

is some small number ot , This number C (cr) is then used as a threshold for 

b 

comparison with the actual value of G after each iteration. If C is 
greater than Cj^ the process is reiterated on the assumption that it is 
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very improbable that this value of C is in fact a minimum. 

Clearly using such a test will reject a certain number of true mini- 
mum calues of C. For this small percentage (cr) of cases ^ the stopping 
criteria is based upon the relative change in state estimates* When the 
relative change in each component of the state vector after an iteration 
is less than some small number, EPS, the successive estimates are con- 
sidered to be equal and the process is terminated. 

On the other hand, passing this first test does not assure that a 
minimum has been reached. For this reason, a second test is prescribed 
which specifies a minimum improvement* Choosing the minimum improvement, 

C^y may be considered purely arbitrary. On the other hand it may be 
helpful to invoke a statistical interpretation to aid in choosing C^^. 

Such an interpretation exists if all of the random sequences are Gaussian. 
Then C(;c^) is related to the likelihood of and CCr^j-CC^^C^"^^) is related 
to the likelihood ratio. The likelihood ratio is a common statistic used 
to test the significance of the difference between two estimates. 

The difference is considered to be significant at the 3 level if the 
probability of occurrence of the observed likelihood ratios is less than 
3 under the assumption that is the true state* The test of statisti- 
cal significance is made by setting a threshold on the likelihood ratio, 
or some function of it. The difference, , is minus twice the 

likelihood ratio and has a chi square distribution with the number of de- 
grees of freedom equal to the number of components in the state, X. C it 
chosen as that value for which the probability of a chi square variable 

less than C is 0 . The test is: if is greater than C re- 

it L 

iterate, otherwise terminate the procedure. The satisfaction of the test 
suggests the interpretation that the last two estimates are not significantly 
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different so no further iteration is carried out. 



The convergence test described may be summarized in terms of three 
quantities involved in the minimization process and three thresholds. 

The three quantities are 1) r, the maximum absolute relative change in 
any component of the state vector from one iteration to the next, 2) C^, 
the current value of the cost and 3)(C^ the change in the cost over 

the last iteration. The corresponding three threshold parameters are 
EPS, C^, C. . The decision rules for terminating or continuing the pro- 
cess are displayed in Figure 1. 




Figure 1. Flow graph of the criteria for iteration termination. 
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As an example^ consider a single stage minimization where the state 
has six components and the measurement has four components. The pro- 
bability that a chi square variable having four degrees of freedom will 
have a value greater than 14.9 is 0,005. This suggests that if C ^ 14.9 

V* 

only one true minimum out of 200 will be rejected by this test. The prob- 
ability that a chi square variable with six degrees of freedom will be 
less than 0.872 is 0.01. If C(X^)-C{X^ =.872 the last two iterations 

are not significantly different at aO.99 confidence level. Finally, for 
those unusual cases where the true minimum is greater than 14.9 the 
iterations are continued until the iteration values are the same within 
the limitations of computer word length. For the GDC 1604 the floating 
operations carry about ten significant digits. It should be considered 
that absolute convergence has been attained when the relative change in 

all of the components of the estimate do not change more than one part in 

9 -9 

10 , This indicates a choice for EPS of 10 

The remarks of this section were directed toward providing some in- 
sight into the choice of the threshold parameters. While the assumptions 
which would make these interpretations rigorous may in most cases be 
lacking, the filter designer must incorporate a convergence test, i.e., he 
must choose a set of parameters. The interpretations discussed above are 
offered as an aid toward choosing an efficient set of parameters. 
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8. Criteria for the Number of Smoothing Stages to be Carried. 

As a result of the way in which the best estimate has been defined, 
the complexity of the algorithm increases with the number of stages over 
which data is available. The estimate of state at the initial time is 
the result of a minimization involving n (number of components in the state 
vector) variables. The estimate of the state at the time of the twelfth 
measurement would be the result of a minimization over 12n variables. 
Clearly, proceeding in this way the computational requirements will ex- 
ceed the capabilities of any computer after a finite number of stages. 

The work of Denham and Pines [5] has focused attention on the dif- 
ference between the case where the measurements are linear in the states 
and the more general nonlinear case. This difference was shown to be the 
result of neglecting higher-order terms in the expansion of the measure- 
ment function. The minimization over many stages may be viewed as a means 
of avoiding this problem since each linearization is only tentative. As 
new data become available, providing more information about the old states, 
the linearization of the measurement functions becomes more accurate. 

When the state is known well enough so that the second-order terms are 
negligible, the linearized measurement is considered to be accurate. This 
approach will be developed into a criterion for the number of smoothing 
stages which must be carried in the next minimization process. 

Consider a second-order expansion of a single nonlinear observation 
function about a point * 

Z = h(x°) + h h Cx°) (jj-of °) + U (91) 

X XX 



where h 0f°) 

X^ and h 
o 

at X • 



is a row vector of partial derivatives of h evaluated at 
) is the symmetric matrix of second partials of h evaluated 
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In order to form a synthetic observation, of the form of (16), 
the synthetic observation must be a linear function of the true state 
with added noise which has zero mean* 

Z° ^ Z ^ h(jf°) + h 0c°) ;e° - b (92) 

X 

= /( je + y + 1>' (93) 

where ^ = h (x^) and 
X 

b = h ErCx-;f°)^h (x°)0f-jf°)] (94) 

XX 

and v' is the variation of the second -order term about its mean^b. 

The added noise D ’ is considered to be the random noise caused by 
the linearizing process* The magnitude of this nonlinear noise can be 
measured in terms of its variance, Ji' . Expressions for evaluating b 
and fi' are developed in Appendix II* 

The process of dropping a stage of smoothing is a lumping operation* 
It can be seen by examining the equations for the linear filter that all 
past data is brought forward in time through (17) and (21)* This is not 
done imnediately in the nonlinear filter because the observation equations 
have been linearized at a point which may be quite different from the true 
state* By keeping several stages active in the filter it is possible to 
perform the linearization at a point much closer to the true state* The 
dropping of a stage should be accompanied by a high degree of confidence 
that the last linearization was performed at a point close to the true 
state* That is, H is not going to change significantly as better esti- 
mates of the true state become available* The invariance of H is related 

to the expression for 7?* = %:race[h Of°)P] since h (x°) represents the 

XX XX 

variability of H with x and P represents the variance of Thus when 

the natural noise in each element of the observation vector for a given 
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stage is an order of magnitude greater than the nonlinear noise, the 
stage corresponding to observation no longer carried. 

The mechanics of lumping are quite straight forward once it has 
been decided that the current linearization is a good approximation. 
Under these conditions, the Kalman sequential filter equations are dir- 
ectly applicable. If these equations are applied to those stages which 
are to be lumped, the result will be an a priori estimate of the first 
state which is not to be lumped. The states which have been lumped no 
longer appear. All of the information that these states carried with 
respect to the estimation of future states is characterized by the a 
priori estimate of the first state which is not lumped. From this point 
then, the problem has exactly the same structure as the problem before 
lumping except that there are fewer stages being carried. 
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9. Computation of the Solution. 

The basic features of the algorithm are shown in Figure 2. The 
following is a brief resume of the quantities which must be computed in 
order to implement the filter. Consider the overshoot decision. In or- 
der to decide whether an overshoot has occurred it will be necessary to 
evaluate the cost, C. For this purpose the multi-stage expression (12) 
is more convenient than the single stage expression (12*). To test for 
convergence it is necessary to have the current cost, the previous cost, 
the current estimate, the previous estimate and the two parameters C^ 
and which are functions of the number of smoothing stages currently 
being carried. The decision on the number of stages to be carried de- 
pends upon an evaluation of the nonlinear noise associated with each 
observation. In order to evaluate this nonlinear noise, the matrix of 
second partial derivatives of the observation functions must be computed 
at the current best estimate of the true state. In addition, there must 
be available a covariance matrix representing the uncertainty of this 
estimate. 



Computation of the Cost 

In general there are many means of computing a given quantity. It 
turns out that the expressions for some quantities are useful in discus- 
sing the problem but are not the most efficient in actual computation. 
Such is the case with the cost C* For discussion purposes the cost was 
expressed in terms of the single-stage state variable. For computing 
the cost at an actual estimate, however, the form of (12) is more effic- 
ient. 

The first term is handled as follows: 

|lje(l/0) - jf(l/k)|lJ^ = ||b0c( 1/0) - je(l/k))||2 (95) 
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Figure 2. Flow graph of overall procedure. 
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where Wj^ = P^i\/0) (the pseudo inverse of P(l/0)) and B is chosen so that 
the equality holds which implies that B satisfies 

B*‘*B = P^(l/0) . (96) 

A suitable B is found by a special routine which begins by decomposing 
P(l/0) into the form 

P(l/0) = (97) 

where A is an n x r matrix and r is the rank of P(l/0)* If r = n then 

P*iUO) = p"^(l/0) and (98) 

B = a"^ . (99) 

If r £ n then 

P*(UQ) = a''^^A^ (100) 

# # # T -1 T 

which implies that B s A and A can be computed by A = (A A) A , where 

T 

the indicated inverse is known to exist by construction* That is, A A is 
r X r and has rank r so it is non-singular* The routine which computes A 
also computes A ^ if it exists, with only minor additional labor* 

For most applications P(l/0) will not be singular even though the 
single-stage covariance will be singular if T has rank less than the 
system order* It is for this reason that the procedure adopted has a 
built-in flexibility to handle the singular case but handles the non- 
singular case with virtually no loss of computational efficiency over the 
more conventional approach of inverting directly the covariance matrix 
P(l/0)* 

Evaluating the typical second term in (12) is accomplished in a 
similar fashion 

lU(i/k) - tr(i-l/k)||J^ = ||s^;c(i/k) - S 2 x:(i-l/k)|l ^ (101) 

T # 

where ^2 = CT F ) • Since T and ^ are assumed to be constant throughout 
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the problem Sj and S 2 can be computed beforehand. By comparing terms 

s^s^»(rr'^)^ (102) 

and $2 = . (103) 

There is no loss of generality in assuming that T has full rank, i.e., 
rank equal to Its minimum dimension. Under these conditions 

(T = T^V (104) 

where = (T^D’ V . This Implies that 

= d^)" V (105) 

and S 2 = (rt')" (106) 

where the indicated inverse is known to exist. 

Finally the third typical term of (12) is 

l| 2 (k)-hCx(i/k))|lJ (107) 

'^3 



where 

The procedure used to evaluate this term assumes that the measurement 

errors are independent, i.e., Ji is diagonal. While this is a realistic 

assumption for most real problems there are means similar to those al~ 

ready used to handle cases where the observation errors at any time are 

correlated with one another, i.e., is not diagonal. The details will 

not be discussed for lack of physical motivation. 

This computation is carried out in the computer subroutine COST. 

The auxiliary matrices B, and $2 are computed outside the subroutine. 

T 

The matrix decomposition routine which generates A such that AA ■ P(l/0) 
is called DECOMA. 
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The Basic Minimization Procedure 



Another procedure which is computed by a different method than that 
used in the development is the minimization process (43). The computa- 
tional procedure is a step-by-step solution of the linearized single- 
stage problem. This single-stage problem is viewed as involving a se- 
quence of observations. Each observation is combined in turn to produce 
a new estimate of the state and a new covariance of that state. Figure 
3 shows the details of the step-by-step procedure. At any point in the 
procedure the best estii^ate, given the data processed, is x and the assoc- 
iated covariance is p. This type of sequential processing is valid under 
the assumption that the observation errors are mutually independent. In 
the computer program this process is carried out with the observations 
divided into blocks according to the time of the observation. The sub- 
routine KALFIL processes each block in the manner indicated by Figure 3. 




Figure 3. Flow graph of step-by-step minimization procedure. 
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This routine is entered once for each block or once for each smoothing 
stage carried. 

This process yields the useful result that the vector x after j 

blocks is composed of the sequences of state estimates jf(l/j), Jr(2/j), 

2 

je(k/j). Similarly the matrix P is composed of k submatrices of the 

form 

P(l,l/j) P(l,2/j) . . P(l,k/j) 

P = P(2,l/j) ♦ ♦ « . ♦ • 

P(k,l/j) . . . . P(k,k/j) 

It is at this point that a lumping operation takes place. If it has been 

c h 

decided that all of the observations up to and including the j stage 
have negligible nonlinear noise, then a lumping operation is performed. 
This consists of shifting all of the estimates x np and out and shifting 
the P matrix up and to the left. This has the effect of eliminating any 
reference to any state at time j or earlier and reducing the dimension of 
the single-stage state x and the single-stage covariance p. 



X = 



'xil/2) 








CM 

CM 

H 


AFTER 


SHIFTING ^ 


x(l/0) 


xO/2) 


SHIFTING 


*(4/2)J TI® INDEX > * - 


x(2/0) 


xi^/2) 









p = 



P(l,l/2) P(l,2/2) P(l,3/2) P(l,4/2) 

P(2,l/2) P(2,2/2) P(2,3/2) P(2,4/2) 

P(3.1/2) P(3,2/2) P(3,3/2) P(3,4/2) 
P(4,l/2) PiU,2/2) P(4,3/2) P(4,4/2) 



AFTE 



SHIFTI 



h-p- 



>(3,3/2) P(3,4/2)| 


SHIFTING 




[pa, 1 / 0 ) 


P( 1,2/0) 


_P(4,3/2) P(4,4/2^ 


TIME INDH 


[p(2,l/0) 


P(2,2/0) 



Figure 4 Evolution of the estimate and covariance during the transition 
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The process of shifting in the computer produces an automatic change in 
indices. The process is shown as two step for an example where j a 2, 
k a 4 in Figure 4. The subroutine SHIFT performs the details of shifting 
the matrices as indicated above as well as several other matrices which 
must be shifted. 
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10. A Simple Example. 



Unfortunately examples seem to fall in two mutually exclusive cate* 
gorles, enlightening or realistic. The following example is introduced 
to illustrate the mechanical details of the algorithm. This example 
also illustrates the process of abstracting the mathematical model from 
the physical situation. Consider an active, drunken, tight rope walker. 

He is put on a tight rope so that only one coordinate will be needed to 
specify his position which will be designated r(k). A. drunken person is 
considered in order to introduce the concept that his position is a ran- 
dom quantity. 

It will be assumed from previous experience with drunken tight rope 
walkers that his next position is different from his last position by 
some completely random variable. The mean squared value of this dif- 
ference is assumed to be known. Further it is assumed that his ramblings 
to the right are balanced in the long run by those to the left, l.e., they 
have zero mean. The mathematical model for this much of the physical 
situation is given below. 



For concreteness Q is taken as 0.02. 

The first expression is usually said to be the model of a dynamic 
system excited by white noise. The second and third are quantitative de- 
scriptions of the noise. 

The next part of the situation that must be described is the process 
by which data are obtained. Here the concept of a nonlinear measurement 
is introduced. An angular measurement is made at some fixed distance 



x(k+l) = x(k) + ou(k) 
E[u)(k)j = 0 



(108) 



(109) 




( 110 ) 
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from the tight rope, i.e*, the observer turns his head through a certain 
angle. For simplicity the observer is placed opposite the center of the 
tight rope at a distance of one unit* It will be assumed that the ob- 
server can sense the angular deflection of his head with a standard de- 
viation of 0.1 radians. The mathematical model for the observer is given 
below. 

^(k) = tan'^(je(k)) + U(k) (111) 

E[u(k)] = 0 (112) 

E[»(k).,(J)J - III 12 (113) 

and B has been taken as 0*01* 

Finally the a priori data must be specified. For this example the 
use of an artificial a prior i will be illustrated. The physical situation 
is such that before taking any data there is essentially no Information 
available about position of the tight rope walker. This fact is model- 
led by taking the a priori estimate as zero and assigning a very large 
variance to this estimate, say 10,000. 

The structure is 

x(l/0) - 0 (114) 

P(l/0) - 10,000 (115) 

This completes the mathematical model of the physical process underlying 
the observations. Assume that the first two observations are 0.7854 and 
0.900 radians. Using these observations the computations required by 
the filter are described in detail* 

The method described in Section 5 for the single-stage case is 

applied to this example for the first observation. The first estimate 

of the state is the a priori estimate “ 0. The partial derivative 

1 12 

of the measurement is evaluated to find H » l/[l+(^j^) ] ■ 1 for the 
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first Iteration. Applying (44) yields 

- (1)-(10,000)/[(1)-(10,000)»(1) + (oOl)] - 1. 

Next is computed from (40) , 

- .7854- tan" ^(0) + - 0.?854o 

From (43) 

Tj (not an exponent) ■ 0 + (1)'(7854 - * 0.7854 

In Table I the results of repeating this procedure three times are shown. 

Also shown in the table is the cost associated with each estimate, 

including the cost for the a priori. From a table for the chi square 

variable the threshold values are found to be C, (0.05) “ 3.84 and 

C (0.95)“ 0.004. Both of these are for a single degree of freedom. C, 

has one degree of freedom since there are two constraints » the a priori 

and the observation; less one adjustable quantity* the single component 

of the state estimate. The also has a single degree of freedom since 

M 

there is only one element in the state vector c From Table I it can be 

2 

seen that the cost associated with ^ might be considered a minimum cost, 
but there has been a significant decrease in the cost so the process is 
repeated. Considering the last iteration it may te said that IcO is not 
a significantly better estimate of the true state than 0«9767« This 
terminates the first minimization process r 

It is interesting to note that the device of taking a large a priori 
variance has led to the expected result that converges rapidly 

to tan [2^(1)] * 1.0, It may occur to the rei^der that this is the hard 
way to evaluate tan [^(1)3, but the advantage of general applicability 
of the method outweighs the advantages of considering special cases. In 
any case, the machinery for handling a priori information must, be avail- 
able in order to implement the lumping procedure. 
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TABLE 1 



EXAMPLE OF SINGLE-STAGE MINIMIZATION PROCESS 



A PRIORI ESTIMATE - 0.0 OBSERVATION - .7854 

A PRIORI VARIANCE - 10,000 OBSERVATION 

BRRC» - .01 
VARIANCE 



ITERA- 

TION 

NUMBER 

1 


OLD 

ESTI- 

MATE 

1 

X 

1 






i 

z 


NEW 

ESTI- 

MATE 

i+1 

"^1 


COST 

C(0?^) 


1 


0.000 


1.000 


1.000 


0.7854 


0.7854 


61.68 


2 


0.7854 


0.6185 


1.617 


0.6042 


0.9767 


1.40 


3 


0.9767 


0.5118 


1.954 


0.5121 


1.0000 


0.015 


4 


1.0000 










.00001 
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Next the bias and variance of the nonlinear noise are evaluated to 
determine whether a lumping operation is indicated. The nonlinear noise 
depends upon the variance of the new estimate and the second derivative 
of the observation function evaluated at the estimate. The variance of 
this estimate is computed from (54): 

» 10,000 - (10,000)^/(40,000 + 0.01) = 0,04. 

Using this variance the lumping test criterion can be computed from 
Appendix II. The bias due to the second order terms is 0.01 and the vari- 
ance of the second order term is 0.0001. Comparing this with the variance 
of natural noise (observation errors) it is noted that there is an order 
of magnitude difference and a lumping operation would normally take place. 

For this example the first stage will not be lumped so that the de- 
tails of the two- stage minimization process may be Illustrated. If 
lumping had taken place the estimate would be projected forward to form 
the a priori estimate for the next time frame. Since 1 for this 

simple dynamic system the new a priori is just the old best estimate. 

From (21) the variance of the a priori estimate is .06. To obtain the 
estimate of the position of the tight rope walker at the second time 
frame one proceeds exactly as above using the new a priori information. 

In order to solve the two- stage minimization problem the problem 
is reduced to a single-stage problem. The elements of the single-stage 
state are the position at the first time and the position at the second 
time. The relations described in Section 6 are used to generate a com- 
plete single-stage problem having a two-dimensional state. 



— - 




— — 


u 




X 




£ 


^o 


0 
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p - 

o 



10000 

10000 



10000 

10000.02 






0.7854 

0.9000 



(117) 



(118) 



P 




0 




(119) 



h(jc) 



hj^(jf) 


a 


-1/ X 

tan (jfj) 


h2(Jf) 




tan ^(^ 2 ) 



( 120 ) 



The above is a new single-stage minimization problem and conceptually, it 
is solved in exactly the same way that the previous (one dimensional) pro- 
blem was solved, l.e., by repeated application of (43). As a practical 
matter it is expedient to adjust the estimate separately for each ele- 
ment in 2. This is possible since the errors in the observations are in- 
dependent. This is true in general because the errors were assumed to be 
white in the multi-stage problem. 

Each Iteration proceeds in two steps. First both elements of the 
estimate are adjusted for the first element of 2 and then the resulting 
adjusted estimates are adjusted for the second element of 2 . See Figure 
3. The result of the first adjustment is already known and need not be 
computed. The first element of this intermediate result is the best 
estimate from the previous minimization process. The second element is 
just the predicted estimate Jf(2/1) = (1/1) ” 1.0 from (83). The in- 

termediate covariance is computed from (85), (86), and (87) and known 
value of P(l/1). _ 



1 

^ inter. 



1.0 

1.0 



( 121 ) 
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Inter. 



0.04 0.04 



0.04 0.06 



( 122 ) 



Starting from this intermediate result the adjustment to the esti- 

1 

mate for the second measurement is computed. H is now a row vector of 
partials of the second observation function with respect to both elements 
of the state vector: -Co. 0.5oJ . Note that the zero in ^ is a 

general result of the way in which the problem is formulated. The measure- 
ment function is formally a function of the whole state x although it is 
clear from the construction of the single-stage form that each individual 
element of the measurement function, h(.r), is a function of the state of 
the system at only one time. From (44) the adjustment coefficients, 
are obtained. 

^1 



.8 

1.2 



(123) 



Applying (40) the synthetic observation is found to be 
= .9 - tan‘^(l) + (.5). (1) « .6146 



and the adjusted estimates are 



1.0917 

1.1375 



(124) 



The cost for the intermediate estimate (121) and the above (124) 

2 

estimate, are 1.313 and .552 respectively. Since these are a sum of 

four squared residuals with two adjustable quantities the convergence— 

test parameters are different. C (.05) ■ 5.99 and C (.05) ■ .103. Based 

L N 

on these parameters it may be inferred that the above estimate is signifi- 
cantly better than the intermediate estimate. A second iteration will be 
carried out. 
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The minimization is accomplished in two steps. The first step is to 
adjust the a priori vector for the first observation. The difference be- 
tween this step and the single stage minimization is that the partial 

2 

derivatives are evaluated at given above in (124). The result of 
this step is comparable to the previous intermediate estimate in (121) 



and (122) 



2 

^ inter. 



.9951 

.9951 



(125) 



p 

inter. 



.048 


.048 


.048 


.068 



(126) 



The second step is to adjust this intermediate estimate for the second 

2 

observation. The partial derivative, ff, is to be evaluated at and not 
2 

X . ^ . The result of this second adjustment is 

inter 



3 

X 

1 



1.0978 

1.1405 



(127) 



The cost evaluated at this estimate is .539. The convergence tests in- 
dicates that the last estimate is not significantly better than the 
previous one. This minimization is said to have converged. 

After it has been decided that the process has converged it is neces- 
sary to evaluate the second-order terms in the expansion of the nonlinear 
measurement function. For this purpose it is necessary to have the co- 
variance of the last estimate. This covariance is automatically computed 
by the method displayed in Figure 3. 



.02096 

.02096 



.02096 

.02966 



(128) 
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From Appendix 11 the expected value of the second order term, denoted 
as the bias or simply b, is 



b 



.0048 

.0065 



9 



and the variance about this expected value is 



(129) 



.000023 



.000042 



(130) 



Assume that only the first of the measurements has negligible second 

order terms. Then a lumping operation is Indicated. This might be ac* 

2 

compllshed by noting the two-stage interpretation of x ^inter* 



inter . 



X(l/1) 

x(2/l) 



(131) 



and 



inter. 



(132) 



P(l,l/1) P(l,2/1) 

P(2,l/1) P(2,2/l^ 

2 ~ 

The lower element of x. ^ and the lower-right element of P, ^ consti- 

Inter ® inter 

tute the a priori Information for the state at the second time frame. 

It would be possible to consider this a priori information and the second 
observation as a new problem. 

In the computer program it is inconvenient to store all of the 
intermediate results awaiting a decision on which stages are to be lumped 
An alternate method for carrying out the lumping operation will be de- 
scribed. Assume as above that it has been decided to lump the first 
stage. The next steps in the filter operation would normally be as fol- 
lows. The third state is predicted using (83) and the covariance matrix 
augmenO>d accordingly using (84). These estimates are adjusted for the 
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third observation. The main minimization procedure Is begun. Recall 
that the main minimization procedure Is carried out In three steps, an 
adjustment for each observation. After the first of these steps the the 
Intermediate result can be Interpreted as 



JC(1/1) 

Jf(2/1) 

xO/l) 

and 

P(l,l/1) P(l,2/1) P(l,3/1) 

P(2,l/1) P(2,2/l) P(2,3/l) . 

P(3,l/1) P(3,2/l) P(3,3/l)^ 

At this point the lumping operation Is carried out by reducing the dimen- 
sion of the single stage to two (see Figure 2) and storing the lower part 

of X. (133) and the lower-right part of P (134) In the area as- 
Inter Inter 

signed to a priori Information. The remaining two steps of the main mini- 
mization procedure are then carried out. If the process has not converged 
then the next minimization will only have two steps. 

This completes the description of the operation of the filter for 
this very simple example. 



Inter. 



(133) 



(134) 



■^Inter . 
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11. Target Tracking. 

It was decided to exercise the scheme on as realistic a problem as 
could be found. The target data (which was fed into the filter) was gen* 
erated by a sophisticated simulator. The target motion is the result of 
maneuver commands generated by the user of the simulation scheme. The 
simulator then computes the motion of the target, and simulates the radar 
returns which that target would generate. The simulated radar is of the 
search type having as available outputs range, range rate and three dir* 
ection cosines at a rate of one frame every two seconds. The simulator 
decides, taking into account the relative position of the target and radar, 
whether a return is received. If a return is received, the simulator out* 
puts a noisy version of the true range, range rate and direction cosines. 

If no return is received the simulator sets a flag in the output data. 

At long ranges the chances of getting a return are relatively small but as 
the range decreases the radar gets returns more and more consistently. 

Forming the Mathematical Model 

It should be noted that this problem, as sketched above, does not fall 
directly into the model which has been assumed in the development of the 
technique. Among the parameters which have not been given in the descrip* 
tion of the problem are r, B and even x (the state space). This is 
typical of the way in which a problem is first encountered. What follows 
will be a series of engineering approximations which yield the mathematical 
model. This model forms the basis for the filter design. 

First consider the stochastic dynamic model. The dynamic model may 
be viewed as specifying two features of the problem. The first is a pre- 
diction function. One asks: how would one predict some future state of 

the system given perfect knowledge of the present state? This question in 
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fact helps to define the concept of state. The state of the system (for 
filtering purposes) is that collection of current attributes of the system 
which has a bearing on the future of the system. For the aircraft target 
the assumption of straight and level flight leads to an assignment of pos- 
ition and velocity as the states. The prediction function is based upon 
the assumption of constant velocity. Thus the components of the state 
vector are 

= north position (miles) 

= north velocity (miles/sec) 

- east position (miles) 

a east velocity (mlles/sec) 

- down position (miles) 

= down velocity (miles/sec) 

and the prediction function is linear in the states and of the form 
jf(k + (h) . $ is the discrete time form of three Independent 

double integrators for a sample time of 2 seconds. 



1 2 



0 0 



0 0 



= 




0 



0 0 



1 2 



0 0 




0 0 



0 0 



0 0 



1 2 



0 0 ' 0 0 ' 0 



1 



The second feature which the model must provide is a measure of the 

prediction errors, or equivalently T* This implies that the prediction 

errors are random variables made up of several normalized random variables. 

T 

The prediction error covariance is then C = F F • For this problem 
F was chosen under the assumption that in each direction there would be a 
step-wise constant component of acceleration of random amplitude having 
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zero mean and variance CT • This implied that F takes the form 




The parameters a^j» and must be chosen to account for turbulence and 
pilot maneuvers, which, from the filter's point of view, appear as random 
accelerations. 



a, = .02 

ajj= .0032 

Secondly consider the observation process. Having chosen the state 

space X it is straightforward to write the functions h(v) 

2 2 2 h 

hj^Cc) = ’^3 ^^5) = range (miles) 

h2(x) = 2 2 2 ~h — “ range rate (miles/sec) 

Ofi + 



0?1 

hoCr) = — 5 ;; FT = north direction cosine 

^ OfJ + 0^3 + xly 



^3 

h^(x) = — 2 2 2 ~^ “ east direction cosine 

CXi + Jf3 + 



X5 

hcCx) = — 2 2 FT “ direction cosine 

Op^ + ;?3 + xp^ 

The noise variance was approximated from considerations related to the 
radar simulator. 
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1.5 



0 



0 



0 



0 



0 



n = 



0 



0 

0 



4.0 X 10 



0 

0 

0 



0 

1.8 X lO’^ 
0 
0 



0 

0 

3.1 X 10"^ 



0 



0 

0 

0 

w A 

1.5 X 10 ^ 



Finally, the a priori estimate and its associated covariance must be 
specified in order to complete the mathematical model. The fact that the 
first radar return has been received places the target in a certain volume 
in physical space. The position components of the a priori estimate were 
taken as the centoid of that volume and the limits of the volume were taken 
as three standard deviations on either side of the centroid. The a priori 
velocity estimates were based on the assumption that the target was headed 
directly toward the radar (in the negative north direction). The magnitude 
of the velocity was taken as that of a Mach 2 target. The covariance of 
each velocity estimate was assumed to be large compared to the square of 
this velocity. The a priori estimate was taken as: 



:e(i/0) = 



80.0 

-0.3 

0 

0 

0 

. 0 



The a priori covariance was taken as: 
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p(i/0) . 



400 

0 

0 

0 



0 

.09 

0 

0 



0 

0 

400 

0 



0 

0 

0 

.09 



0 

0 

0 

0 



0 

0 

0 

0 



0 0 
0 0 



0 0 400 0 

0 0 0 .09 



This completes the process of abstracting the physical situation into 
the form of the mathematical model. 

It should be emphasized that the abstraction process must be carried 
out for each physical system that generates a sequence of measurements. If 
the model accurately describes the conditions under which the measurements 
are made then the filter can be expected to yield estimates which are best 
in some sense. Even an accurate model and an optlmimi filter do not assure 
that the estimates will be adequate for any particular purpose. 



The Algorithm Parameters 

There are three parameters that define the iteration termination cri- 
teria. They are the probability, a, that a minimum cost is greater than a 
given threshold, C^; the level of statistical significance, 3; and the num- 
ber of significant digits used by the computer. 

The threshold, C, , depends upon the number of stages carried (which 
determines the ntimber of degrees of freedom) as well as upon a. There are 
five degrees of freedom in the cost for each stage carried since the system 
dynamics (1) Introduce six constraints and the observations (5) introduce 
five constraints and there are but six adjustable parameters (the components 
of the state vector) for each stage carried. An a of 0.05 was chosen in 
order to have a small but finite number of cases where the minimum cost 
was greater than C^. 
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The likelihood- ratio- test threshold, C^, also depends upon the num- 
her of stages carried. The number of degrees of freedom is six for each 
stage carried since that is the number of components in the state vector. 
A significance level of 0.95 was arbitrarily chosen. The thresholds C 

Li 

and C are shown in Table II. 

M 

The filter was implemented on a CDC 1604 computer. This machine 
carries about ten significant figures. Two successive iterations were 
considered to be equivalent if all of the components of the estimate were 
equal in the first nine significant figures. 

The lumping criterion is based on a comparison of the covariance of 
the observation errors and the covariance of the nonlinear noise intro- 
duced by the linearization process. For concreteness, the nonlinear 
noise was considered to be negligible when its covariance was less than 
that of the observation errors by a factor of ten. 

Target Tracking Results 

Three target trajectories were filtered. The results were quite 
similar. The target on which the largest number of observations were re- 
ceived will be described in some detail. 

Figures 5 through 9 are a graphical display of the filter operation. 
Figures 5 and 6 show the true target trajectory projected on the NORTH- 
EAST plane and the NORTH-DOWN plane. Superimposed on the true trajectory 
are confidence areas generated by the filter. The boxes are used to pro- 
vide a measure of the quality of the estimate. The size and shape of the 
box is computed from the covariance matrix of the estimate. If the es- 
timation errors were Gaussian with a covariance equal to that computed 
by the filter the box would have the following interpretation. There is 
an ellipse, centered about the estimate which contains the true state 
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TABLE II 



ITERATION CONTROL PARAMETERS AS A FUNCTION OF THE NUMBER 
OF STAGES CARRIED FOR O' «.U5 AND 3 = .95 



NUMBER OF 
STAGES CARRIED 


^L 




1 


11.07 


1.64 


2 


18.31 


5.23 


3 


25.00 


9.39 


4 


31.41 


13.85 


5 


37.65 


18.49 


6 


43.77 


23.02 


7 


49.55 


27.86 


8 


55.76 


32.85 


9 


61.33 


37.80 



68 




[I3] 




Fig. 5. Projection on NORTH-EAST 
plane of true trajectory 
with estimates. 



Fig. 6. Projection on NORTH-DOWN 
plane of true trajectory 
with estimates. 
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c 



Fig. 7. Observation On) and estimation (+) errors in 

NORTH coordinate. 
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Fig. 8. Observation (]□) and estimation (+) errors in 

EAST coordinate. 
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with probability of 0.63 and whose boundary is a curve of constant prob- 
ability density. The vertices of the box are the extremities of the 
major and minor axes of that ellipse. 

Figures 7, 8, and 9 show the estimation errors and the observation 
errors in each of the position coordinates. In addition the computed 
standard deviation of the estimates is shown as a solid curve. 

In order to illustrate the ability of the algorithm to converge to 

a leas t-^squares estimate the cost is given in Table III as a function of 

the number of iterations and the number of observations received. The 

first cost listed in each row was evaluated at an estimate which has 

not been adjusted for the newly received observation. The estimate 

for a (k+l)-stage minimization process is given by (89). The first iter- 

2 1 

ation yields x^ by adjusting X^ for the most recently received observation. 
Only the partial derivatives associated with this new observation are 
evaluated for this step. This step is comparable to the simple single- 
stage linearization employed in [4] with such disappointing results. The 
second row indicates the danger of stopping at this point. Subsequent 
iterations reevaluate all of the partial derivatives at the previous 
estimate. All of the observations are reprocessed with estimates start- 
ing from the a priori estimate. 

The cost after a lumping operation is only the sum of squares of 
the residuals associated with stages still carried by the filter. The 
lumping operation occurs in the middle of second iteration because this 
is the first time that the intermediate results needed to form the new 
a priori estimate of the remaining stage become available again after the 
decision to lump has been made. 

The time required to produce a leas twsquares estimate is tabulated 
in Table IV. The time indicated does not include the time required for 



72 



TABLE III 



TYPICAL SEQUENCE OF COSTS AS A FUNCTION OF NUMBER 
ITERATIONS AND NUMBER OF OBSERVATIONS PROCESSED 



OBSERVATION 

NUMBER 


C(Jfl) 


STAGES 

CARRIED 


1-1 


1-2 


II 


1®4 


1 


367. 


0.57 


0.13 




1 


2 


377.8 


118.5 


2.07 


2.07 


2 


3 


35.4 


11.7 


11.16 




3 


4 


19.0 


12.9 


12.9 




4 


5 


48.7 


20.99 


20.99 




5 


6 


121.6 


27.36 


27.36 




6 


7 


811.9 


38.28 


38.06 




7 


8 


749.9 


101.3 


44.16 


44.16 


8 


9 


48.51 


44.4 


9.49* 


9.49 


2 


10 


12.27 


10.98 


1.51* 


1.51 


1 




* A lumping operation occurred. 
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TABLE IV 



TYPICAL COMPUTATIONAL TIME REQUIREMENTS FOR 
THE CDC 1604 COMPUTER 



OBSER- 

VATION 

NUMBER 


COMPU- 

TATION 

TIME 

(SEC) 


CUMULATIVE 

COMPUTATION 

TIME 

(SEC) 


OBSERVATION 

ARRIVAL 

TIME 

(SIMULATED) 


1 


0.817 


0.817 


0.0 


2 


2.633 


36.633 


34.00 


3 


2.617 


48.617 


46.00 


4 


3.683 


52.300 


48.00 


5 


5.633 


65.633 


60.00 


6 


7.633 


73.266 


62.00 


7 


10.45 


83.716 


66.00 


8 


26.083 


109.800 


78.00 


9 


14.900 


114.700 


80.00 


10 


2.100 


116.800 


82.00 


11 


1.367 


118.167 


86.00 


22 


1.467 


134.050 


110.00 


34 


1.300 


148.417 


130.00 


63 


1.483 


196.017 


196.00 



74 



auxiliary computations performed for diagnostic purposes nor the time 
required to read the data from the magnetic tape* It does include all 
computations inherent in the filter operation such as evaluating the 
cost and covariance of the nonlinear noise. The cumulative running time 
has been adjusted to reflect the fact that the filter cannot begin to 
compute a new estimate until the next observation is available. 

Comparison between Table III and Table IV yields the obvious fact 
that the time required to compute the estimate is highly dependent on 
the number of stages carried. From an analysis of the computations in- 
volved in Figure 3 it can be shown that the computations increase as the 
square of the number of stages times the system order. The computation 
time depends linearly on the number of observations. 

The operation of the filter indicated for the tenth observation is 
typical of all of the remaining stages in both number of iterations and 
processing time required with the exception of the cases where the minimum 
cost was greater than C. . This happened 3 times out of 67 observations on 
the longest run. In each case additional iterations involved only a single 
stage. Two or three extra iterations were needed to satisfy the termina- 
tion criterion. The largest relative change in any component of the es- 
timate decreased approximately two orders of magnitude after each itera- 
tion. The average processing time for these three cases was about 1.9 
seconds . 
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12. Conclusions 



An algorithm has been developed for the processing of a sequence of 
noisy, nonlinear observations made on a dynamic system whose state is a 
random function of time. The best estimate of the state of the system 
at each observation time is defined to be the weighted-least-squares es- 
timate • 

This estimate is computed by solving a sequence of linear problems 
which approximate the nonlinear problem more and more closely. A method 
has been developed for automatically determining the number of iterations 
required to compute the least-squares estimate by the above procedure. 

The computation of each estimate is based on a least-squares fit on 
only a finite sequence of past observations. A method has been developed 
for determining the length of this sequence of past observations. The in- 
formation contained in the older observations is carried forward in the 
form of an a priori estimate. 

The radar tracking problem is an example of the type of problem which 
falls within the scope of this investigation. The algorithm was implemented 
on a digital computer and used to process a sequence of observations pro- 
vided by a realistic radar-target simulator. The estimation errors were 
generally within the expected range, considering the randomness of the 
dynamic system and the observation errors. The algorithm achieved the 
least-squares estimate in three or four iterations. The length of the 
sequence of observations on which the least-squares fit was based, rapidly 
settled to only a single previous observation. The computational require- 
ments appear excessive when compared with those associated with linear 
observations. There are, however, no other generally applicable methods 
when the observations are nonlinear and there is little prior information 
about the state of the system. 
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Appendix I 



Discussion of singularity of 

The covariance matrix P is assumed to be non*singuIar throughout 

o 

the development, however^ there is a meaningful interpretation for the 

case where P is singular. It will be shown that the filtering equations 
o 

are still valid in view of this interpretation and that in the one 

instance where the inverse of P is required (in evaluating C) the use of 

o 

the pseudo inverse is appropriate. 

This development begins by defining a new set of state variables » 
so that the errors in the a priori estimates are uncorrelated. 

^ = ite (1) 

yo = '^o 

-1 T 

where U is a unitary matrix such that U = U and 

UP = D (3) 

o 

where D is a diagonal matrix. 

It will be shown now that D is the covariance of the a priori 

estimates in the y states . 

Gov iy^) 5 E[ 

Cov (y ) " E[U0f";f ) 

o o o 

Cov (y^) = UE[ Oc’xJ 

Cov (y^) =s D (4) 

If p is singular then D has at least one zero on the diagonal or, 
o 

to be more specific, the rank of P is equal to the rank of D which is 

o 

the number of non-zero elements along the diagonal of D. There is no 
loss in generality in assuming that all of the non-zero elements of D are 
in the upper part of the diagonal. The upper elements of y are 
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conventional statistical estimates of the corresponding elements of the 

true state having a variance given by the element in D. On the other 

hand, the lower elements of are precise or exact estimates of the 

corresponding true state components and have no variance. When these 

interpretations are reflected back to the x states the meaning of a 

singular P becomes clear. A singular P implies that a certain number 
o o 

of linear combinations of the states are known exactly. 

It will now be shown that the filtering equations, used repeatedly 
in the minimization process, produce estimates consistent with these 
interpretations. That is, the estimate of those linear combinations of 
the states which were known exactly before adjustment for observed data 
are not affected by the adjustment process. Further the covariance 
matrix, P^ of the adjusted estimates reflects the fact that these linear 
combinations are known exactly. This demonstration will be carried out 
using the y state space. The filter equations are transformed to operate 
on the y coordinates. 

Consider first 



PJi\hph'^ + [z - Hx^l , 

substituting y for x in (5) yields 

+ PP'^rHPH'^ + Iz - , 

pre-multiplying both sides of (6) yields 

substituting U*^DU for P^ from (3) 

y^= [z - , 

_ T 

and finally making a change of variables f( = H\3 yields the filtering 
equation, 



(5) 



( 6 ) 



(7) 



( 8 ) 



(9) 
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In the y states. Since the lower elements of D are zero It Is clear that 
there Is no change In these components of the adjusted estimate, y^ , 
after application of (9). 

Now consider the covariance equation 

?1 = Po “ ^ ^ 

* U*^DU - + P]”^Pu’^DU 

p^ = u''^[D - + p]’Vd]u (lo) 

T 

multiplying In front by U and In back by U 

up^u’’^ = D - {KVK + p]"Vd (ii) 

Thus P^ reflects the fact that those linear combinations of x which 
were known exactly are still known exactly. 

In the expression for the cost consider only the first term, 

2 

IIjc ~ JiTiIL, * where W. was assumed to be the Inverse of P If It existed. 

" o l"Wj^ 1 o 

Substituting In the y states yields 

•J* i i 

Defining Q = UWU and substituting above yields 

■ -lllw, - 11^0 - 

Define y = ~ y^) and partition “y Into two subvectors 




where the lower subvector y^= 0 from (9) 

jL 
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Let Q be partitioned as 



Q - 




where is the upper, non-zero, part of D. Thus sum of the residuals 
can be expressed as 





I 



The blank submatrices in Q above are immaterial since they will be 
multiplied by = 0. Arbitrarily assigning zero submatrices to the 
blanks implies that those residuals which are known to be zero are given 
zero weight* This leads to 




T 

and preraultiplying by U and postmultiplying by U yields 




0 



0 



U 



but this is just an expression for the pseudo inverse of p . 

o 

Thus if the pseudo inverse of is used in the definition of the 
cost there will be no weighting of the residuals in certain linear 
combinations of the states* On the other hand, if the is always 
computed using (5) or one of its derivatives then these particular 
residuals will always be zero* 

During the discussion of the minimization algorithm the non* 
singularity was an important assumption. The discussion is valid for 
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the case of singular In Che sense that a new minimization problem can 

= . The discussion Chen Implies 

o *1 

Chat Che change In Che estimate of has a positive component In Che 
direction of the negative gradient of C with respect y^. When these 
conclusions are reflected back Co Che x state space It can be seen Chat 
Che change In Che estimate Is related Co the projection of Che gradient 
of C Into Chat subspace of the state space about which there Is some 
uncertainty and for which an adjustment In Che estimate Is meaningful. 



be defined In terms of y^. Caking 
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Appendix II 



Consider a random variable 

0 ■ ( 1 ) 

where A is a symmetric matrix and ;r is a Gaussian random variable 

with mean zero and covariance P . It is desired to find an expression 

for the first and second moment of the random variable v. The 

development begins by making a change of variables 

X • BUj^ (2) 

T 

where B is a decomposition of P such that P « BB , U is a unitary matrix 

T 

as yet unspecified, such that UU = I and ^ is a random vector of zero 
mean and identity covariance. The random variable u is expressed in 
terms of y by substituting (2) in (1). 

y « ^/’^^B’^ABUi/ (3) 

Now let U be chosen so that 

U^b\bU ■ D (4) 

where D is a diagonal matrix. Substituting (4) in (3) yields 

u > y^Dy (5) 

or 'O can be expressed in terms of the components of y 

V - IAv\ ( 6 ) 

i 1 ^ 

To evaluate the first moment of v the order of summation and 
expectation is interchanged. 

E[y] = E[E d yh 
i ^ 

E[U] - 
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The components of y all have unit variance* 

Eru] « ^ d. o) 

i 

Expressing (7) in matrix form yields 

E[u] = trD (7') 

2 

The second moment is evaluated by expanding V in terms of the 
components of j/. 

= Ef (r d d z/^)] 

i j J J 

Efy^] = Efy 
Er.2] , -r 

The first term is evaluated by recalling that the fourth central 

moment of a unit -variance Gaussian random variable is 3. Each 

element of the second term can be factored due to the independence 

of the components of y. 

E[u^] = 31 dj + E d.d.E[2/J] E[y^j 
i i i J i J 

Efy^] = 3T df + r d.d. 

i ^ ^ J 

Efu^] = 2T dj + (T d 
i i 

Substituting (7) above yields 

E[y^] = 2E dj + (Efy])^ (9) 

i ^ 

This implies that the variance of u is 

Var[y] = d^ (10) 

i ^ 

or in matrix form 
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Var[u] = 2tr[D^] . (10*) 

The results, (7) and (10), are expressed in terms of the parameters 
of the original problem. Substituting (4) in (7') yields 

Efu] = tr[U^B^A B U ] . 

Taking advantage of the fact that trfAB] = tr[BA] yields 

E[y] = tr[A B U U^B^] 



and cancelling the unitary matrices yields 

Efu] » tr[AP] . 

The development for the variance precedes along analogous lines* 
Var[u] = 2trf U^B^A B U U^B^A B U ] 

Varfu] = tr[A B B^A B B^] 



( 11 ) 



Var[u] = tr[APAP] 



( 12 ) 



It is possible to consider the covariance of two random variables 
as follows 

V = x"^A X (1) 



and a second random variable 

v' ^ x^k'x (!') 

The means of both of these random variables are known from (11) and an 
expession will be developed for the expected value of the product. After 
making the same change of variables as before, the expected value of the 
product can be expressed as 



E[u'u] = E[/c y /d y1 


(13) 


where 




T T 

C a U B A'B U 


(14) 


This is equivalent to 




E[u'u] = tr(C E\y y^D y y^J) 


(15) 
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Examining only the elements of the matrix Inside the expected value 

T 

operator, it is noted that the middle term, y Ti y% is a scalar. 

y^TiyZ d y\ (16) 

i ^ ^ 

Thus , the elements of the matrix are 

This expected value is zero for i^j and for l’*j it can be vnritten as 



E[j/J E \yh “ E[d + E[i/^ E d 1^] 

IrKK 11 

s d^;^] - 3dj + £^d^ 
k k 

Thus, the expected value Is a diagonal matrix of the following form 



(18) 



E[i/J/^Dl/l7] » 2D + (E d )•! 

k 



Substituting in (15) yields 



e[u v } - tr C[2D + (S d )i] 

k ^ 



which can further be simplified to 



e[u u] - 2tr [CD] + [E d JtrC 

k 

Considering only the factors of the last term, it is noted that 



(19) 



Ed, - e[u] (21) 

k 

and that 

tr C - tr[U^B^A B u] 
tr C “ tr[A P3 

tr C - E[y ] (22) 

Now the first term can be identified with the covariance from the general 
expression 
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(23) 



Cov ir^x^'] - E[x^'} 

Cov[u'y] - 2tr[c d] 

CovCu'y] - 2tr[U^B^A'B U d] 

Cov[y'y] ■ 2tr[A*B U D U^B^] 

Cov [y'y] “ 2tr[A'B U u"^B^A B U uV] 

Cov[y'y] = 2tr[A'P A P] (24) 

The useful results of this analysis are (11) and (24) since (12) is 
only a special case of (24). Since these results will be used as a 
guide even when the random variable x is not known to be Gaussian it is 
worth reviewing just where the assumptions were used. The factor 2 
which appears in (24) comes directly from Gaussian assumption (i*e., the 
fourth central moment is 3 times the second central moment squared). A 
second result of the Gaussian assumption is that the new random variables 
y are statistically independent (used in (8) and (18))l The variables 
y are uncorrelated (since D is diagonal) but only in the Gaussian case 
does this imply independence. The effect of this assumption is diffi- 
cult to evaluate. On the other hand (8) does not depend upon the 
Gaussian assumption. 
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Appendix III 



Program Listings 
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COMMON PHI(8»8)»DEL(8*8)*Q{8.8)*X10(8»10)»P10{80,80),R(8)»S10(8»8) 
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X(l*2) IS USED TO STORE THE TOTAL ( CUMULAT I VE ) COST OVER ALL STAGES 
MAIN FILTER LOOP STARTS HERE 



DO 30 KM=1»MAX 
DO 33 I=1»NS 
33 RR( I*K)=R( I ) 

READ TAPE 3»ITMAX,JT 

READ TAPE 3* GX ♦ ( A , I =2 ♦ I TMAX ) »GY ♦ ( A ♦ I =2 ♦ I TMAX ) ,GZ 
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51 X2 ( I »K)=X2 ( I .K+1 ) 

GO TO 30 
20 CONTINUE 

IGNORE FRAMES WITH NEGATIVE RANGE 
IF (Z(1»K) )21»21.2A 
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AG=AG* 5.-X(l,2) 

AG =(CUBERTF(X( 1»1 )/AG)-l.+ 2 . / ( AG*9 . ) ) /SQRTF ( 2 . / ( 9 « *AG ) ) 

COMPUTE the NORMAL APPROX I MAT l OM TO THE CHI SQUARE VARIABLE 
DISTIBUTI0N» THIS IS USED TO EVALUATE THE OVERALL CONSISTANCY 
OF THE RECIEVED DATA WITH THE MODEL ASSUMED 



OUTPUT INFORMATION ON FILTERED ESTIMATES FOR GRAPHICAL DISPLAY 
ON FRAMES WITH DATA ONLY 
STORE INFORMATION FOR FILTER SUMMARY 
WRITE TAPE 49» NFIL ,KM»K ,SSQAV»AG 
WRITE TAPE 49» ( Z ( I ♦ K ) ♦ I = 1 , 5 ) 
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121 FORMAT ( 32X1 3HSTATE SUMMAR Y ♦ / / , 1 1 X5HNOR T H ♦ 1 2X AHEAST ♦ 1 IXAHDOWN » 
19X10HVELOCITIES»9XAHTIME»//7X32HFILT REAL OBSR FILT RFAL OBSR 
2»32HFILT REAL FILT NORTH EAST DOWN »//) 

302 READ TAPE 49 
READ TAPE 49 
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READ TAPE A9,GX»GY*G 
RF=XF**2+YF**2+ZF**2 
RF=SQRTF( RF ) 
DNF=XF/RF 
DEF=YF/RF 
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PRINT 150 ♦ TIMA ♦ TIN 
307 PRINT 125*KT,KS»AER»COSTA 
125 FORMAT! 5X ♦ 2 I 6 , 2F 1 0 . 4 ♦ / ) 
END 
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THIS SUBROUTINE HANDLES THE FIRST ITERATION ON EACH NEW PIECE OF 
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this subroutine handles the MFCHANICS of the main iteration loop on 
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JU = JUMP( JU ) 

GO TO (9 tlAtZAtZO) , JU 
9 DO 11 J=1»K 
DO 11 I=1»NS 

11 X2( I ♦J)=(X1 ( I »J)+X2( I » J) )* 



CALL AUG MENT ( PIO, XIO) 
CALL AUG MENT ( PB,X2 ) 
K=K+1 

14 DO 15 1 = 1 ,K 



c 

vO 



-- o 

f-H 

•^3 X 

3 “) 





































CM 




















O 














o 


X 




















f-H 






O 








f-H 


•> 




















a 






CNJ 








CL 


CD 






z 






— 








II 
















CL 






o 






»— » 








— ■ 






o 










w 






»— » 






•> 








1 — 1 






CO 


z 














1— 




•> 


3 
















QC 






1— 


1— 


cn 




< 


CO 


3 






0- 


a. 


3 






(\J 


3 






z 


Z 






M 


z 


w» 


O 


CM 


X 


X 


w- 






f-H 


K 




f\J 


LU 


LU 


(\j 




f-H 




(\J 


f-H 


l- 






CO 








UJ 




1— 


z 


z 






51 


f-H 


X 


X 


co 


f-H 




CL 


f-H 


CO 


CM 


QC 




CO 








LU 


1 


II 


II 


II 


O 


II 


II 


II 


1 


z 


f-H 


3 


LU 


O 


o 


o 


1- 


3 


Z 


3 




—X 


u 


»— 1 


3 


—X 


1— 


3 


X— 


^ 3 


3 


U 


3 


3 


f— H 


Z 


1— f 




f-H 


1— • 


II 






3 


f— H 


oc 




Z 


Z 


II 


< 


< Q 


z 


»— 1 


z 


m 




•> 


f-H 


v£» 


v£> 




z 




o 


f-H f“* 


►—» 


f-H 




Z 


w- 


H- 




f-H 


3 


3 


1— 


f-H 


f-H 


►— 1 


II 


3 


1— 


1 H- 


1“ 


1— 


3 


3 LU 




z 


z 






X— 


CO 






X-- 


1— 


3 




X z 


z 


CO 


3 


3 


u. 


o 


f— * 


O •-H 


CM 


o 


o 


o 


CO 


1— • 


< 


o 


II o 


o 


O 


< 


< 


*-• 


u < 


Q 


X 


X 


u 


a 


Q 


a 


z u 


o 


X o 


u 


u 


U 


U 



f\j m vo o CO o 

f-H i-H Cn (N 



u 



98 






cn 

2 

D 

or 



UJ 



=) 

o 

q: 

CO 

3 

to 



cr 


















C 


CO C 
















U- CO 


f—i 3 
















►-H 


ao ^ ^ 
















CO 


^ CO e 
















Z 1- 


O -- X 
















C X 


rH r— 1 <t 
















►— « ►-*< 


to fsl X 
















to I 


0>. ^ ^ 
















*— 1 CO 


X 
















> 


QC OC Z 
















o o 


• — ^ 
















cr u j 


cr 00 3: 
















a cr 


^ z: 










— 






►— « 


^ to 
















to 3 


O -J CO 










Cf't 






Z O 


CO Lu Z 










• 






»— t LL' 


•> o ^ 










o 






< cr 


o 










f-H 






1“ 


00 ^ o 










UJ 






Z to 


^ o 00 
















o 


O rH 


















f-H •> O 


















Q. 00 00 
















e z 


^ w 
















z — 


-- 00 X 










o 






< a 


o a 










u. 






X 


*— < X 




— 






a: 






Z 3 






o 






cr 






O -J 


CO r- O 




i-H 






3 






1— i 


. o 




— ' 






U 






1 - z 


o ^ 




X 






u 






< UJ 


•-H 00 




CO 






o 






cr X 


X 00 -- 




z 












U jS. 


w CM 










z 






h- 


— MX 




— 






c 






►— • 


CO 




c 












to 


— CM 




f-H 






H- 






U H- 


CO 00 O H- 










< 






— z 


^ f-H CO 




cr 






cr 






CO UJ 


o cr o 




z 






UJ 






< 21 


CO 00 LJ 










a 






CD UJ 


«-<. »> w 9^ 




— 






o 






QC 


00 I-H I-H 




o 












UJ 


00 X K- 




f-H 






o 






X 3 


00 to 




9> 






z 






O 


^ oo -- O 




00 






f-H 






UJ 


o u 




#> 






a 






CO cr 


UJ o ^ 




00 






X 








cz t/) X 




— * 




O 


3 






o 


#> 00 <( 




o 




II 


-J 




— 


LU Z 


— — — X 




u 




U 






CO 


z 


00 O < 




CD 




»— 1 


X 




►— 1 


0. 


^ rH NJ ♦> 










o 




— ■ 


^ X 


CO 3 




— 






CO CO 




-J 


3 3 


-- 00 — _l 




o 




3 


X 




h—t 


O -J 


o 


— 


f-H 


1— 


-J 


• — ^ 


to 


Lu 


cr 


X X f-H h“ 


o 


9> 


o 


II 


f-H 


-J 




CO UJ 


CL O 


rH 


00 


o 


CO 


M 


★ 


-J 


3 X 


— 00 O 


#> 




X 


-J 


f—i 


u 


< 


CO 1“ 


Z 00 X 


00 


00 


CO 


(A 


1- 


►— f 


X 




o CO to 






-J 




<T f-H 


1 




to o 


X O -- -J 


cr 


u 




X 


X CM 


»— • 


-J 


z cz 


X ^ 


cr 


CD 


H 


II 


cr 


II 


-J 


X UJ 


O >- X X 




< 


3 


to o o 


CD 


< 


1- CO -J 


u 




9^ 


-J 


X 


X o 


►— » 


u 


z ^ 


I-H (M CO 




lT» 












UJ < 












(M 




m 


to u 












f-H 




CM 



r\J 

cn 

f\j 

(\j 



►-I ll) 
I r> 
D z: 



U- O 
u 

CO 

rg 



CO 

LU 



CO 

UJ 



< 

CO 



O 



X 



II 

=) 



vO 

<\l 

o 

o 



1“ -- Z) 

U Z) H- -J 
Z -J CO ^ 

3 O K 
Li- U to 

rH O f\J 

q: X o u ^ 

^ Lu i-H 

UJ .-H a + 

Z rsj 5 : I- 

♦ 3 -- Z 

^ f\J -J f— •— 

z X ct: 

O UJ *-H h“ Q_ 

Z X h- -- X 

I- < X UJ 

UJ < -J II X to 3 

I- X 3 -- CO Z Z 

<C X )<c i-i 

3-J-J 

-J -J 3 .-H -J II II Z 

<<u--<uao 

>UUXU»-^!yCU 
UJ < 

nO 

CM CM 



o 



u u u 



u u 



99 



21 CONTINUE 
LU= 

LSHOOT = 
END 



this routine initializes the atate vector according to the statement 
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114 FORMAT ( 30H RANDOM VECTOR FROM N(0»I) HAT) 

115 FORMAT ( 30H INITIAL VALUE OF STATE VECTOR) 

116 EORMAT ( 20H SQUARE ROOT OE Q ) 

117 FORMAT { 20H SQUARE ROOT OF R ) 
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THIS ROUTINE READS IN AND SUMMARIZES THE PROBLEM STATEMENT DATA 
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PRINT no 

READ 104, ( R( I ) , 1=1 ,NM) 
PRINT 104,(R(I),I=1,NM) 
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cr 


c 


IK 
















































rH 


z> 










z 




































oc 




z 










3 






































00 












+ 




































c 




X 










3 




































r— ; 


rH 


< 










ii 




































to 


rsl 












< 
































< 


















3 
































X 




— 














kPt 
































h- 




a 


oc 


z 










to 






















U- 














w 


»> 












z 






















h- 










</) 




ex: 


oc 


3: 
































< 














0<‘ 




Z 










r-H 






















z 










UJ 




— 


to 


IK 










ii 






















1—4 










-J 




c 


-J 


to 






UJ 




3 






















h- 














oc 


U-! 


z 






h- 




<!■ 






















to 










u? 






o 








< 




f-H 






















UJ 










1— 1 




o 


» 


— 






1— 














UJ 


























00 


— 


o 






to 




o 










to 












u. 










UJ 






o 


00 










o 










*-4 












c 










if) 




c 


f-H 








-J 




kf) 










0 








vO 














>— 1 


UJ 


rH 


•> 


o 






< 




z 










z 








f-H 




ex 










c 


to 


ex 


00 


00 






3 




3 
















z 


#>. 




UJ 










z 


»— 1 






w 






O 


X 


+ 










ex 






0 






CD 












o 


— ■ 


00 


CD 






»— t 


»— 1 


1—4 










UJ 


UJ 






vO 




Z 










cc 


z 


c 




CL 






> 


ex 


ii 










0 


u 




1— 


rH 




3 










< 




rH 


X 






— 


1— 1 


h- 


< 










ex 


z 




< 






Z 










UJ 


-J 










o 


Q 


< 


1—4 






3 




0 


< 


to 


ex 


0 














z 


< 


00 


— 


o 




f-H 


Z 


5: 








» 






1—4 


UJ 


UJ 


CsJ 




UJ 










►—I 


ex 




o 


r-H 




W- 


»— 4 




to 






(\i 




0 


ex 


h- 


0 


— 




X 










-J 


X) 


o 


f-H 










UJ 


z 






X 




z 


< 


< 


♦—4 


«— 




h- 










z 


h- 


i-H 




00 




if) 


z 


u 












0 


ex 


z 


to 


3 














o 


< 


X 


00 






z 


< 


z 


f-H 






00 




u 


< 


1—4 


z 






UJ 










z 


z 






OsJ 




•> 




< 


il 










UJ 


> 


i— 


0 


3 




X 


















X 




r— 


Ll 


1— I 


1—4 






00 




to 


0 


to 


u 


w 




H 










UJ 




or 








C CsJ 


O 


ex 














u 


U_| 




ex 














X 






r— 


^ (Ni 




rH f-H 




< 








to 




u. 






ex 


ex 




0 










h- 




00 


00 


O h- 




w wr 


X 


> 


f-H 






z 




0 


UJ 


u. 


UJ 


1 




h- 


















f-H if) 




ex o 


1—4 


o 














i— 


0 3 
















U- 


UJ 


o 


ex 


O 




z > 


ex 


u 


o 






2: 




X 


< 




z 


3 




3 










l-H 


X 




to 


00 U 




#> 


h- 




o 


— 




z 




»— 4 


z 


z 


3 






< 












1— 


— 










< 


3 


to 


< 








ex 


»— 4 


0 




3 




3 










UJ 




oc 


— 


f-H f-H 




O Ovi 




3 


to 


3 




h- 




h- 


h- 


1—4 


UJ 


w 




0 










UJ 


LL 




CO 


X 1- 




f-H f-H 




< 


z 






CL 




< 


to 


i— 


h- 


CsJ 




UJ 










if) 


o 


00 




•> to 




^ w 


-- UJ ex 


1 


< 




^ to 




z 


3 


u < 


ex 




















00 


— o 




00 O 


00 U 


UJ 


to 


»— 4 




(\j < 








UJ 


z 


★ 




UJ 










o 


•> 


-J 


-W 


o u 




X 


z 


> 


z 






ex 




UJ 


3 


3 


1—4 


z 




3 










I— 


z 


UJ 


o 


f-H ^ 




00 ^ 


00 < 


o 


★ 


CD 




^ OD 


u. 


u 


< 


3 




u 




3 


h- 










u 


O to 


X 




<w — . 


»— 4 




3 


a 




CD 


0 


z 


3 


0 to 


w 




< 


tO 








h- 






•> 


00 < 




O (>sJ 


(\i ex 


UJ 


ii 


il 




^ UJ 




< 


0 


u 


UJ 


U- 




> 


UJ 








if) 




— 


— 


— 51 




U rH 


ex < 


X 


z 


r— 




h- 


X 


K- 4 


♦—4 






1—4 






1— 


rH 




rH 


UJ 


o 


00 o 


< 




CD 


> 


h- 


3 


1—4 




ex < 


1—4 


ex 


> 


3 


UJ 


to 




< 




CsJ 




CsJ 


h- 




r. 


rH 


fSl 




h- 


— o 




to 






UJ 2 : q: 


< 


*— 1 


3 


X 








UJ 






•> 




1— 


00 




3 




^ •— 1 


00 U 


iL 




3 




0 


h- 


> 


0 


0 1- 


z 




Z 


X CsJ 




0 


UJ 


< 




00 


-- -J 




O 




o 




W 




ex 1 - 


< 


0 


Z 


X 




z 




0 1- 


(\J 




CsJ 


z 


ex 


1—4 




O 


*— 


rH Z 


00 UJ ex 




h- 


f-H 


0 to 


Z 


u 


»— 4 


:* to 














» 






X X 


f-H h- 


O 


o 


— X 


u. 


f-H 


CL 


1 


UJ 










UJ 


rH 


i— 


to 


to rH 




rH 


h- 


z 


CL 


#> 


o 


f-H 


00 *-• 


I— H* 




ii 


il 


3 


u 


UJ 


UJ 


UJ 


UJ 


h- 


ii 


to 


UJ 


LO 


CsJ 




CsJ 


3 


•—1 






00 O 




if) 


CL 


UJ 


3 




ii 


UJ UJ 


X 


X 


X 


X 


< 


3 UJ 


X 


< 






— 


o 


< 


z 


00 


X 


00 


00 Z 


• h” 


h- 




3 




to X h- 


i— 


i— 


h- 


u 




h- 


< 


CL 


tH 


z 


00 


(X 


h" 


o 


w 


00 to 


<w 


UJ 


— o u 


< 


in 




CL 


H- 










r-4 


0 




h- 




1 


ex: 


1 




ex 


5: o 


-J 


ex 


U Z 


00 rH < 


z 


f-H 


*— 4 


2: 


3 


to 


to 


to 


to 


Q 


Cvi U 




X 


3 


3 




to 


UJ 


z 


f-H 


rH ^ 


ex 


CD 


-- ti ex 


1—4 




»— 


3 


3 to 


1—4 


»— 4 


1—4 


»— 4 


z 




1—4 


CL 


u 


<— 


1— 




*— « 


u 


o 


>- 


X ^ 


•> 


< O 


CD Z h“ 


h- 


O 


h- 


3 


< 1-4 










1—4 


0 to 


z 


1—4 


U- 


UJ 


u. 


X 




u 










U X 


to 


o 


CL 




U 


Ovi 


Ovi 


h- 


CsJ 




0 < 


3 


X 


»— 4 


ex 


1—4 


1— 


< 




f-H 


Ovi CO 




m 


f-H UJ 


UJ 








CD ex 


ex 


CL 


X 


3 




CD 


3 






















































0 


f-H 


CsJ 






















f-H 
























rH 


CsJ 


CsJ 


u 


u 












u 


u 








u 


U 


u 


U 


u 


u 




U 


U 


u 









107 



20 CONTINUE 

SET UP TRANSFER OF NONLINEAR TERMS TO BASIC FILTER 



DO 18 J=1,NM SZA( J,L) =R2( J,J) 
18 Z1 ( J*L)=Z1( J,L)-8( J) 

15 CONTINUE 
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1 .F10,2»10H REL CHANG ,F10.7) 

2 JUMP = 2$PRINT 102,COST2*A $RET(jpN 
102 FORMAT(40H REPROCESS COST2= 

1 ♦F10.2»10H REL CHANG ♦F10.7) 
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This routine evaluates the sum of the squares of the residuals 
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0 6 L=1,NS 
= A+XA(L»K)* PHI ( I jL) 
XA(I,K+1) =A 
END 
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14 BCD( J3»J4*J1 )=BCD( J3»J4»J1+LU) 

X10( J4»J1 )=X2( J4,J1 ) =X2 ( J4, Jl+LU) 

4 X1(J4,J1)=X1(J4,J1+LU) 

DO 15 J1 = 1,K $ NSK ( J1 ) =NSK ( Jl + L'J) 

15 NR( J1 ) =NR( Jl+LU) 
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+X( 3* 



A(5»5,1)=1./R+X(5»K)*H3 
A(1,3,1)=A(3,1,1)=H1*X(3,K) 
A(1,5,1)=A(5»1»1)=H5*X(1»K) 
A( 5,3, 1 )=A( 3,5, 1 )=H3*X ( 5,K) 






ro 














,-H 


CO 


CO 


lA 


lA 


w 








X 

1 


X 

1 


X 

1 


X 

1 


X 

1 


X 

1 


< 

II 








1 

</) 


1 


I 

CD 


CD 


1 

CD 


1 

CD 








q: 


(T 


q: 


or 


or 


or 


<r 












* 


* 














(/) 


lA 


CD 


CD 


CD 


CD 


,-H 








(X 


QC 


QC 


ct: 


IX 


or 














★ 




★ 




lA 








CO 


lA 


1 


lA 


rH 


CO 


-w 








X 


X 


X 


X 


X 


X 


< 








★ 






★ 






II 








r-H 


,-H 




CO 


lA 


lA 










X 


X 


X 


X 


X 


X 


CO 








★ 


★ 


★ 


★ 


★ 


★ 












r-H 


CO 


CO 


lA 


lA 


lA 








X 


X 


X 


X 


X 


X 














* 


★ 


* 




(O 


• 


• 


• 


• 


• 


• 


• 


• 


• 




CO 


CO 


CO 


CO 


CO 


CO 


CO 


CO 


CO 


< 




★ 


* 


II 


II 


II 


II 


M 


II 


II 










--*> 












r— H 


CO 


lA 


CO 


CO» 






lA 


lA 


lA 



X X X 
I I \ t—i 1—i 
Ll Ll Ll 
ct ct: q: CO m 



CO CO lA lT\ 



lA 



CO CO 



cococo<<< <<<< 



u 

♦—I 

CD 



★ 


★ 




II 


II 


II 


II 


II II 


II 


















< 






★ 


* 


★ 






— 


— • 






CD 












1- 




*— • 






— - 


— >* 


— 


CO 


CO 






lA lA 




X 












X 




























★ 












★ 




< 












CO 


lA 


rH 


lA 


rH CO 


lA 


CD 
















•-H 






r-H 


CO 


lA 














X 














u 


— 






«—> 




«— 


rH 


rH 


CO 


CO 


lA lA 


rH 


★ 




CD 


CD 




CD 




II 


rH 






X 


K 


X 








— 


^ — 


w 


CO 21 




Z 


z 




z 


u 




X 


X 






«— 




< 


< 


< 


< 


< < 


< 


X Z 






•N 








*— • 


+ 


z 


*— i 


II 


1? 


II 


II 


II 


II 


II 


II II 


II 


★ 




rH 


rH 




rH 


#« 










— ' 










— 


— 






lA rH 




II 


II 




II 


< 


CD 


*—* 


rH 


rH 


CO 




lA 




lA 


CO 


lA 


CO <f 


lA 


X II 


• 


< 


CD 




u 








n 


II 




















)fC *-H 


o 


t—* 


*—* 




*—* 


w 




CD 


*— • 




,-H 


CO 


lA 


rH 


rH 


CO 


CO 


lA CA 


CO 


rH 


II 










< 


< 


II 
























X Csj 




c\} 


CO 


• 




+ 


»— • 


— 


lA 


lA 


rH 


CO 


lA 


rH 


rH 


CO 


CO 


lA lA 


rH 


★ 


I— I 






o 




u 


w 


t—t 






w 






w. 


— 


W 


— 


•w 


— 


• O 




O O 


II 


o 


II 


rH 


*— 


O 


O 


< 


< < 


< 


< < 


< 


< < 


< 


CO O CD 


O 


o 


u 


o 


u 


X 


CD 


o 


o 



CO C>sJ 



vO lA 



118 



DO 6 IA=I ,NS 
DO 6 IB=1,NS 

C=C+R1( IA,IB,I)*R1(IB,IA,J) 

R2( I ,J) =R2( J, I ) =C 

END 



SUBROUTINE SKIP DATA 
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35 X(K,J)=X(K»J)-RATIO*X(L,J) 000330 

36 CONTINUE 0003A0 

34 CONTINUE 000350 

40 DO 43 I=1*N 000360 

I1=N+1-I 000370 
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SUBROUTINE DECOM A ( OA ,8 ♦ C ,NS , NR ♦ EPS »N1 ♦ N2 ) 
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0(N(K)»N(J))=Q(N(J)»N(K))=Q(N(K)*N(J))-B(N(J)»n*B(N(K.),I) 
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B(N( J) , I ) =0(N( J) *N ( I ) )/A 
DO 4 J=IA»NS $ DO 4 K=J»NS 

Q(N(K)*N(J))=0(N(J)*N(K))=0(N(K)»N(J))-B(N(J),I)*B(N(K),I) 
NR = NS 
END 



ENTRY TIME 



CM 



UJ 



or 



c 

c 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o cn 

O 

o r- — • r- 

* O O O t/) O 

^ O W II II II o It 



r- 



“^OU-J'-OO-- 

jz<u>i-2:> 

mUJ-)i/)QLOUJQ 



LU 

2: 



+ + 



rH CM 

f- o 

i-H y 

i/) lO ' 
II II 



(\i cn 



<0-J000~)0 

i-i-<oc:oc:Qi_i2: 



125 



FORMAT( 26 HO INITIAL TIME IN SECONDS ♦ FIO 



o 

f\J 



-- o 

r-H 

-t U_ 



c 

vO 






to 


U- 


















0 


1— 


















Z- 


< 


















0 


0 


















u 


-j 


















LU 


u_ 


















U) 


•f 


















z 




















*— • 


CM 




















1 — 


















UJ 


»— I 


















u 




















z 


U_ 


















UJ 


1 — 


















0:1 


< 


















UJ 


0 


















u_ 


-J 














0 




u. 


u_ 














H- 






+ 














U 




0 


• 














UJ 






0 










C 




0 




UJ 


vO 










h- 








2: ^ 


* 


i—H 








u 








#— • h- 












UJ 


0 






1 — »-l 




CM 








Ci 


h- 














0 




1 


u 






0 U- 




i-H 




1 — 


H- 


1 — 


UJ 






X 1 - 




— * 




u 


U 


U 0 UJ 


H* 


00 < 




• 


• 


UJ 


UJ 


UJ 




CL 


U 


(\J 0 




r- 


0 


a 


0 


a 


* 


< 


UJ 


w _J 






II 




II 


II 


ir\ 


h- 


Q 


h- U- 




1 


• 












11 


< II 




X 


r- 0 


H- 


0 0 


0 


H“ 


UJ 


0 


X H- 




X 


m h- 


Z 


1 - H“ 


h- 


z 


1 — 


1 — 


(X u 




w 


II u 


♦— 1 


U 


u 


»— 1 


♦— « 


U 


0 UJ 




Ll 


X UJ 


a 


UJ 0 


UJ 


a: 


QC 


UJ 


u_ 0 




*— 1 


X 0 


a 


0 0 


a 


Ol 


3 


0 




i-H 


















in 






rH 






CM 









z: 



2: 

o 

u 



o 

z 

UJ 



126 



IDENT RNDEV63 



c c 



: C 


c 




c 


c 


e e 




c 




c 


c 


e 


C 


c 


e 


c 


c 


c c 


c o 


o 


C 




c 


c 


C 


c 


o 


1 cn 




IT 


vO 


r- 


00 c^ 


C 


1— 1 


r\j 


cn 




iT 


o 


h- 


00 


c^ 


c 


I-H (\J 


CO ^ 


ir. 


vO 


1^ 


00 


c^ 


C 


r-H 


CM 








c 




c c 




f—i 


r— 1 


f—i 


I— ( 


f-H 


r—i 


f-H 


i-H 


1—H 


<\J 


(\J CM 


CM CM 


(M 


CM 


CM 


CM 


CM 


CO 




CO 








c 


o 


c o 


c 


o o 


C 


c 


c 


c 


c 


c 


c 


c 


c c 


C C 


c 


O 


c 


c 


c 


c 


c 


c 


^ c 


c 


c 


c 


o 


o o o 


c 


c 


o 


o c 


c 


c 


c 


c 


o 


c c 


c c 


c 


c 


c 


c 


c 


c 


c- 


c 


' c 


c 


c 


c 


c 


o c 


c 


c c 


c 


c o 


c 


c 


c 


c 


c 


c c 


o c 


c 


c 


c 


c 


c 


c 


c 


c 


' c 


c 


c 


c 


o 


o o 


c 


c o 


c 


c c 


c 


o 


o 


c 


c 


o c 


c o 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


o c 


c 


c c 


c 


c c 


c 


c 


c 


c 


c 


c 


o 


c 


c c 


o c 


c 


c 


c 


c 


c 


c 


c 


c 






















































z 






















































o 
















































Ll 






» — • 
















































H- 






H- 
















































< 






< 
























LS 




Q 


C/ 
























*— 1 
























CL 




CL 


or 






X 












> 






> 
























< 




CO 


< 






H- 












LU 






LU 




































vO 












Q 






o 












Lj» Ll 












h- 






h- 






<t 






























*— 1 *— • 












3 




ilJ 


3 


















LU 






d 












z z 












CL 




X 


a 






LU 


KD 










> 






CL 








CL 




3 3 












Z 




»- 


z 






X 


LU 










1— 1 






< 












z z 












1—1 






H- 1 






»- 


z 
















3 








c 




















O 




















< 






Z 








CL 




<D <D 












Ll- 




H- 


Ll 






o 


Lu 










KD 






< 








< 




CL CL 












O 






c 






h- 












LU 






K- 












< < 
















rvj 




















z 






tO 


X 






o 
















to 




Q 


to 






(M 




















LU 






z 




h- K 












to 




C 


to 








z 










CL 






>- 


o 






o 




to to 




»- 








LU 




X 


Ll' 






KD 


LU 










O 






CD 


z 






u 




CL CL 




»— 1 








QC 






CL 






U < 


X 


< 








Ll 








*— 1 






LU 




*-H *— • 




X 








Q 




UJ 


O 






< 3 


LU 
















LU 








to 




Ll Ll 




LU 








3 




u 


C 






CL Ll 


3 


H- 








iw' 






O 


LU 






















< 




3 


< 






H- 


0- 


< 








KJ 






1 — • 


> 






h“ 




H- K 




H- 








H 




a 


II 






CO h- 


X 


o 








Ll» 






> 


< 






LiJ 




LU LU 




LU 








★ 




LU 








3 LU 


o 


3 








X 






f-H 


to 






to 




to to 




to 








★ 




CL 








to to 


u 


Ll 








KJ 






o 



cn cn 

vO vO 

> > 

LU H- LiJ 

O •- O 

Z X z 

(i: uj q: 



»— < 4* 

I CL Q_ 



X ^ O O X 

c UJ (\J QC Q: .-H LU 









I-H 






















CO CO 


CM 








+ 






















4- 4- 


4- 








to 








to 














tO to 


tO 






CO 


z 








Z 














z z 


z 


LU 




i-H 


o 




LU 


£L 


O LU 




LU 


CD 






LU 


o o 


o 


to 




LU 


u 




to 




KJ to 




to 


r- 






to 


u u 


u 


< 




> 


> 




< 


> 


> < 


I-H 


< 


m 






< 


> > 


> 


CL 


to 


I-H 


3 




CL 3 


3 CL 


+ 


CL 




r- 


vO 


CL 


3 3 


3 


LU 


r-H )je 


Lu 


CL 


4* 


LU 


CL 


CL LU 


♦ 


LU 


CM 


<t c 


CO 


LU 


CL CL 


CL 



CsJ 



q: 

H-“>33<—l</)3_|*-i_j<<i-,<.-H_j<QaCQ<Q.U<t0<iO5k'5:2:3 



CO 




vO 




> 


CL 


LU 


U 


3 


> 


Z 


3 


CL 


CL 



127 



STA ** ** = ADDRESS OF SECOND ARC R 00000330 

EXIT ENI 1 ** 00000340 

SLJ 00000350 

ERASE BSS 1 00000360 

FIVE13 DEC 1220703125 00000370 



4* 



c c c 
oc a c 
r' 

c 

c c c c 
c c c c 
c c c o 
c c c c 



c 


o 




c 


o 


lO 


c 


c 


o 


o 


o 


rg 


c 


o 


r- 


c 


c 


fSJ 


c 


o 


r- 


o 


c 


rH 


c 


c 


cn 


c 


c 


fO 


c 


c 


r- 


c 


c 


o 


c 


o 




c 


o 




o 




r- 


(\J 




r*H 



»- h- H- H- O 

u u u u z: 

O O O O LU 



</) 

z 

o 

u 

> 

o 

a: 



128 



INITIAL DISTRIBUTION LIST 



No. Copies 



1. Defense DocumenCaClon Center 20 

Cameron Station 
Alexandria, Virginia 22314 

2. Library 2 

Naval Postgraduate School, Monterey, California 

3. Ordnance System Command HQTS , 1 

Department of the Navy 
Washington, D. C. 20360 

4. Prof Harold A. Titus 1 

Department of Electrical Engineering 

Naval Postgraduate School, Monterey, California 

3. LT Ralph E. Hudson 1 

FAIRECONRON ONE 

FPO San Francisco, California 96667 



129 



--UMC^ASSIFIEB 

Security Classification 



DOCUMENT CONTROL DATA • R&D 

(Smcurity ctmmsHicmtion of tlttm, body of mb 9 tract and Indexing annotation muet be entered whan tha ovarall report ie claeeifled) 



1. ORIGINATING ACTIVITY (Corporate author) 

U. S. Naval Postgraduate School 
Monterey. California 



2a. REPORT SECURITY CLASSIFICATION 

Unclassified 



Zb GROUP 



3. REPORT TITLE 



FILTERING NONLINEAR MEASUREMENTS 



4 DESCRIPTIVE NOTES (Typ» ol nporl and tnclualva dataa) 

Ph D Dissertation. December 1966 



5 AUTHOf^CSJ (Laet name, first name. Initial) 



HUDSON, Ralph E., LT USN 



6 REPO RT DATE 

December 1966 



7a. TOTAL NO. OF PAGES 



130 



7b. NO. OF REFS 



11 . 



6a. CONTRACT OR GRANT NO. 



b. PROJEC T NO. 



9a. ORIGINATOR’S REPORT NUMBER^S; 



9 b. OTHER REPORT NOfS^ (A ny Other numbers diat may be assigned 
this report) 



10 AVAILABILITY/LIMITATION NOTICES 



Qualified requesters may obtain copies of this report from DDC 



n. SUPPLEMENTARY NOTES 



12. SPONSORING MILITARY ACTIVITY 



13 ABSTRACT 

An algorithm is developed for estimating the state of a linear dynamic 
system excited by a random sequence. The input data are noisy observations 
which are nonlinear functions of the state. The estimates are best in the 
sense of least squared residuals. A significant problem in radar tracking 
is investigated and the effectiveness of the algorithm verified. 



DD 1473 



UNCLASSIFIED 



131 



Security Classification 



imCLASSIflED 

Security Classification 



14. 

KEY WORDS 


LINK A 


LINK B 


LINK C 


ROLE 


WT 


ROLE 


WT 


ROLE 


WT 


Filter 
Nonlinear 
Observation 
Radar Tracking 
Least Squares 





























INSTRUCTIONS 



1. ORIGINATING ACTIVITY: Enter the name and address 
of the contractor, subcontractor, grantee. Department of De- 
fense activity or other organization (corporate author) issuing 
the report. 

2a. REPORT SECURITY CLASSIFICATION: Enter the over- 
all security classification of the report. Indicate whether 
'^Restricted Data’* is included. Marking is to be in accord- 
ance with appropriate security regulations. 

2b, GROUP: Automatic downgrading is specified in DoD Di- 
rective 5200. 10 and Armed Forces Industrial Manual. Enter 
the group number. Also, when applicable, show that optional 
markings have been used for Group 3 snd Group 4 as author- 
ized. 

3. REPORT TITLE: Enter the complete report title in all 
capital letters. Titles in all cases should be unclassified. 

If a meaningful title cannot be selected without classifica- 
tion, show title classification in all capitals in parenthesis 
immediately following the title. 

4. DESCRIPTIVE NOTES: If appropriate, enter the type of 
report, e.g., interim, progress, summary, annual, or final. 

Give the inclusive dates when a specific reporting period is 
covered. 

5. AUTHOR(S): Enter the name(s) of authoKs) as shown on 
or in the report. Entcn- last name, first name, middle initial. 

If military, show rank and branch of service. The name of 
the principal author is an absolute minimum requirement. 

6. REPORT DATE: Enter the date of the report as day, 
month, year, or month, year. If more than one date appears 
on the report, use date of publication. 

7a. TOTAL NUMBER OF PAGES: The total page count 
should follow normal pagination procedures, i.e., enter the 
number of pages containing information. 

7b, NUMBER OF REFERENCES: Enter the total number of 
references cited in the report. 

8a. CONTRACT OR GRANT NUMBER: If appropriate, enter 
the applicable number of the contract or grant under which 
the report was written. 

8b, 8c, & 6d, PROJECT NUMBER: Enter the appropriate 
military department identification, such aa project number, 
subproject number, system numbers, task number, etc. 

9a. ORIGINATOR’S REPORT NUMBER(S): Enter the offi- 
cial report number by which the document will be identified 
and controlled by the originating activity. This number muat 
be unique to this report. 

9b. OTHER REPORT NUMBER(S): If the report has been 
assigned any other report numbers (e/ther by the originator 
or by the aponaor), also enter this number(s). 

10. AVAILABILITY/LIMITATION NOTICES: Enter any lim- 
itations on further dissemination of the report, other than those 



imposed by security classification, using standard statements 
such as: 

(1) "Qualified requesters may obtain copies of this 
report from DDC.’’ 

(2) "Foreign announcement and dissemination of this 
report by DDC is not authorized. ’’ 

(3) "U. S. Government agencies may obtain copies of 
this report directly from DDC. Other qualified DDC 
users shall request through 

' . ” 

(4) "U. S. military agencies may obtain copies of this 
report directly from DDC Other qualified users 
shall request through 

9t 



(5) "All distribution of this report is controlled. Qual- 
ified DDC users shall request through 

If the report has been furnished to the Office of Technical 
Services, Department of Commerce, for sale to the public, indi- 
cate this fact and enter the price, if known. 

IL SUPPLEMENTARY NOTES: Use for additional explana- 
tory notes. 

12. SPONSORING MILITARY ACTIVITY: Enter the name of 
the departmental project office or laboratory sponsoring (pa^ 
ing for) the research and development. Include address. 

13. ABSTRACT: Enter an abstract giving a brief and factual 
summary of the document indicative of the report, even though 
it may also appear elsewhere in the body of the technical re- 
port. If additions! space is required, a continuation sheet shall 
be attached. 

It is highly desirable that the abstract of classified reports 
be unclassified. Each paragraph of the abstract shall end with 
an indication of the military security classification of the in- 
formation in the paragraph, represented as (TS), (S). (C), or (U). 

There is no limitation on the length of the abstract. How- 
ever, the suggested length is from 150 to 225 words. 

14. KEY WORDS: Key words are technically meaningful terms 

or short phrases that characterize a report and may be used aa 
index entries for cataloging the report. Key words must be 
selected so that no security classification is required. Identi- 
fiers, such as equipment model designation, trade name, military 
project code name, geographic location, may be used aa key 
words but will be followed by an indication of technical con- 
test. Tbe aasigniment of Unka, and wal^ia la optlxmal. 



DD 1473 (BACK) 



UHCLtSSIFIED 

132 Security Cla«sification 



1 




1 




p 




it 

*1 

\ 



c 



i 



